CN112013792A - 一种复杂大构件机器人面扫描三维重建方法 - Google Patents

一种复杂大构件机器人面扫描三维重建方法 Download PDF

Info

Publication number
CN112013792A
CN112013792A CN202011114453.9A CN202011114453A CN112013792A CN 112013792 A CN112013792 A CN 112013792A CN 202011114453 A CN202011114453 A CN 202011114453A CN 112013792 A CN112013792 A CN 112013792A
Authority
CN
China
Prior art keywords
camera
point
point cloud
dimensional
calibration
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
CN202011114453.9A
Other languages
English (en)
Other versions
CN112013792B (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.)
Nanjing Nanxuan Heya Technology Co ltd
Original Assignee
Nanjing Zhipu Photoelectric Technology 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 Nanjing Zhipu Photoelectric Technology Co ltd filed Critical Nanjing Zhipu Photoelectric Technology Co ltd
Priority to CN202011114453.9A priority Critical patent/CN112013792B/zh
Publication of CN112013792A publication Critical patent/CN112013792A/zh
Application granted granted Critical
Publication of CN112013792B publication Critical patent/CN112013792B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01BMEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00Measuring arrangements characterised by the use of optical techniques
    • G01B11/24Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
    • G01B11/25Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01BMEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00Measuring arrangements characterised by the use of optical techniques
    • G01B11/24Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
    • G01B11/25Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
    • G01B11/2504Calibration devices
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01BMEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B11/00Measuring arrangements characterised by the use of optical techniques
    • G01B11/24Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures
    • G01B11/25Measuring arrangements characterised by the use of optical techniques for measuring contours or curvatures by projecting a pattern, e.g. one or more lines, moiré fringes on the object
    • G01B11/2518Projection by scanning of the object
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00Geometric image transformations in the plane of the image
    • G06T3/40Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4038Image mosaicing, e.g. composing plane images from plane sub-images
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/30Determination of transform parameters for the alignment of images, i.e. image registration
    • G06T7/33Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
    • 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
    • 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/10028Range image; Depth image; 3D point clouds

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明公开了一种复杂大构件机器人面扫描三维重建方法,属于三维形貌测量技术领域。其首先采用格雷码编码图像结合相移的结构光编码方法进行面阵三维形貌测量;接着将面阵三维形貌测量系统和工业机器人结合起来,进行系统标定;最后对采集的点云数据进行处理,进行粗拼接和精拼接处理。本发明的复杂大构件机器人面扫描三维重建方法能够实现高速度和高精度的大型物体三维形貌测量。

Description

一种复杂大构件机器人面扫描三维重建方法
技术领域
本发明涉及复杂大构件机器人面扫描三维重建方法,属于三维形貌测量技术领域。
背景技术
基于面结构光与工业机器人的大构件三维测量,是指把基于编码面结构光的三维测量方法(面阵三维测量)和工业机器人相结合,从而实现大构件物体的三维形貌测量。由于单纯的面阵三维测量方法只能测量固定区域内单个面上的三维形貌信息,通过引入工业机器人,实现一种三维拼接技术,使得面阵三维测量范围不仅局限于单次测量幅面的面积,而且能够测量一个大型结构件的完整三维形貌。
现有的技术在测量过程中存在相位误差,并且在计算过程时,尤其在搜索对应点的过程中,计算量非常大;在寻找对应点集时,可能会产生一定数量的错误对应点,影响最终的匹配结果。
因此,需要一种复杂大构件机器人面扫描三维重建方法。
发明内容
为了解决上述技术问题,本发明公开了一种复杂大构件机器人面扫描三维重建方法,其具体的技术方案包括:
步骤1:采用格雷码结合相移的结构光编码方法进行面阵三维形貌测量;
步骤2:将面阵三维测量系统和工业机器人结合起来,实现系统标定;
步骤3:对采集的点云数据进行处理,进行粗拼接和精拼接处理;
所述步骤1包括依序的双面相机标定、光栅投影与图像同步采集、编码图像的相位展开、基于相位的双目立体匹配和空间点的三维坐标求解;
所述步骤2系统标定为双目标定和手眼标定结合;
所述步骤3包括粗拼接和精拼接两部分,所述粗拼接为使用手眼标定方法,把机器人工作区域内的所有单位三维测量数据统一到机器人基底坐标系下;所述精拼接指的是通过点云处理的方法,把粗拼接存在的微小误差消除掉。
进一步的,所述双面相机标定包括根据相机模型标定出两个相机的内参和外参,确定双目相机
各自的相机坐标系和世界坐标系的组成参数;标定出两个相机坐标系之间的转换关系,得到相应的参数,确定双目视觉测量模型。
所述光栅投影与图像同步采集包括将编码面结构光图案通过投影仪投射到物体表面;投影仪投影时,将正弦条纹图案和格雷码图案按照顺序依次投射出去;此时相机与投影仪同步,每当一幅编码图案投出去后,相机立刻拍摄被物体表面调制的编码图像;然后把图像保存下来,用于后续处理。
所述编码图像的相位展开包括对保存的图像进行解码,获取每幅图像像素点的相位值。
所述基于相位的双目立体匹配包括在编码图像经过解相位之后,左右图像得到每个像素点的相位值;同时根据双目立体视觉原理,遍历图像上的所有像素点,求出左右图像上的所有极线对;然后在极线对的两条极线上,搜索相位值相同的点,即获得左右图像上的匹配点对;因此,根据相位相等条件和极线约束条件,能够实现双目相机里图像像素点的匹配。
所述空间点的三维坐标求解包括根据物点、像点、两个相机相对位置关系以及双目标定的结果参数,在图像二维坐标恢复物体三维坐标下,通过最小二乘法求解物点的三维坐标值。
进一步的,格雷码是一种二进制循环码,其特点是在一组数的编码中,任意相邻的编码值只有一位二进制数不同,并且最大数和最小数之间的编码也只有一位数字不同,因此格雷码又称为循环码或反射码。由于投射的是黑白相间的二进制编码图案,且相邻二进制编码图案之间只有一位二进制不同,使得格雷码编码图案具有极高的可靠性和鲁棒性,也使解码工作十分方便。以7位格雷码为例,需要投射7幅图案,把待测区域分成27个单独的区域,投射的7位格雷码编码图案如图1所示。
格雷码法的基本思想是构造一组具有不同排列方式的二值边缘图案,使投影模式空间可以划分成若干个截面,每个截面都可以通过其在图形序列中的二值强度识别出来,这样就可以对每一个阶段的包裹相位都做上相应的标记从而实现相位展开。基于格雷码的相位展开图如图2所示。
在做格雷码条纹设计时,假设投影仪内所加载的标准正弦相移条纹的分辨率为F, 周期为T,那么所需要设计灰度编码条纹的个数为
Figure 461926DEST_PATH_IMAGE001
,此外格雷码条纹的周 期与相移条纹周期也需要保持一致。
通过将格雷码编码图像加载到投影仪内,经投影仪投射后被相机捕获,条纹解码的过程则通过将所采集的灰度图像经灰度映射外加设定阈值完成。最后,如果某一点像素解码为B(x,y),则该点的绝对相位值为:
Figure 831596DEST_PATH_IMAGE002
(1)
在解码过程中,格雷码编码图像黑白交界处在图像二值化过程可能会存在误差,造成该部分相位展开会出错,出现相位突变情况,在图像中显现毛刺现象,如图3(a)所示。为了获得准确的解包裹后的相位信息,本节主要通过采用自适应中值滤波器来消除解包裹相位误差。中值滤波的模板大小通常选择9×9或者11×11。假设经相移算法获得的包裹相位的分布函数为φ 1 (x,y),采用格雷码直接解码获得绝对相位的分布函数为Φ 1 (x,y),经中值滤波后所获取的绝对相位分布为Φ 2 (x,y),则理想绝对相位Φ(x,y)可通过如下公式获得:
Figure 549016DEST_PATH_IMAGE003
(2)
公式(2)中,round()为四舍五入函数,最终所获得理想相位便是解码准确时所获得绝对相位,消除误差后的绝对相位分布情况如图3(b)所示。
进一步的,所述双目标定和手眼标定结合包括:
首先:硬件设定:设置棋盘格相对于机器人的底座固定不动,且始终位于双目相机的公共视场内,相机捕获棋盘格图像;
其次:双目标定:控制工业机器人带动固定在末端的双目相机,以不同的位置姿态拍摄棋盘格的标定图片;双目相机标定获得到各自相机的内参、外参、两个相机之间的转换矩阵和相机1到相机2的转换矩阵记为M
最后:手眼标定,在进行手眼标定时,把双目相机看成两个没有联系的相机,再分别对各自单独的相机进行手眼标定;棋盘格手眼标定方法是针对单目相机,得出两个单独相机的各自手眼标定参数,相机1的手眼转换矩阵记为X 1 ,相机2的手眼转换矩阵记为X 2 ;对于一个相机1的相机坐标系下的三维点P 1 ,其在相机2的相机坐标系下的三维点P 2 和在机器人末端坐标系下的三维点P e ,有如下两组转换关系:
Figure 334439DEST_PATH_IMAGE004
(3)
Figure 352073DEST_PATH_IMAGE005
(4)
公式(3)中,先把P 1 点转到相机2的坐标系下,然后再相机2的手眼标定矩阵转换到机器人末端坐标系下;在公式(4)中,直接把P 1 点通过相机1的手眼标定矩阵转换到机器人末端坐标系下;由上面两个等式可得:
Figure 95907DEST_PATH_IMAGE006
(5)
在理想情况下,相机1和相机2的手眼标定矩阵存在如上关系;但在实际情况下,由于存在机器人绝对定位误差等原因,使得上式不完全相等,即存在手眼标定误差;因此,可以通过上式的约束条件来对手眼标定矩阵进行误差校正;在双目面阵三维测量过程中,一般把三维测量数据所在的三维坐标系定义为相机1的相机坐标系;由于把相机1的手眼标定结果作为三维拼接的转换矩阵,用相机2的手眼标定结果对相机1的手眼标定结果进行误差矫正;通过以下计算方式,能够得到更加准确的相机1手眼标定矩阵,如公式(6)所示:
Figure 300623DEST_PATH_IMAGE007
(6)
其中,X 0 即为通过相机2矫正后的相机1手眼标定矩阵。
进一步的,测量系统是基于工业机器人,由于工业机器人普遍存在绝对定位误差等问题,使得单纯的通过手眼标定方法实现三维拼接存在一定的误差。为实现高精度的大构件三维测量,本节的三维拼接方法分为粗拼接和精拼接两部分进行。粗拼接指的是使用手眼标定方法,把机器人工作区域内的所有单位三维测量数据统一到机器人基底坐标系下。精拼接指的是通过点云处理的方法,把粗拼接存在的微小误差消除掉。
对于基于点云配准的精拼接方法,使用的是基于点云滤波的ICP点云配准方法,具体包括两个部分,即点云滤波和点云配准,具体流程如图4所示。在点云配准之前,需要进行点云滤波等预处理。由于两个单位点云数据之间的粗拼接误差比较小,所以点云配准方法采用比较经典的ICP点云配准算法。下面分别点云滤波方法和ICP点云配准算法分别进行介绍。
ICP点云配准算法虽然精度很高,但是也有相应缺点,即耗时很长和对匹配初值敏感。为了提高ICP算法的效率和稳定性,有必要在配准前对点云进行预处理。通过使用点云滤波的方法,来改善ICP点云配准的效果。
点云滤波,和图像滤波、信号滤波等相似,是指滤除点云数据中不必要的部分,为后续的点云处理提供便利。点云滤波方法主要有直通滤波器、体素格滤波器、统计滤波器、半径滤波器等,通常组合使用完成任务。下面为各个点云滤波器的定义:
(1)点云直通滤波器
点云直通滤波器,是指在点云的XYZ轴的某个方向上,通过设置固定的数值区间,删除数值区间内的点云数据或者删除数值区间外的点云数据,从而达到保留某些区域点云数据的目的。直通滤波器一般用于点云截取等过程。
(2)点云体素格滤波器
点云体素格滤波器,是指在点云空间内用大小一样的体素格(相当于小方块)填满,每个体素格内最多取一个点,从而达到点云均匀降采样的目的。通过设置体素格的不同尺寸大小,可以实现不同点云数量的降采样。体素格滤波器对一般用于稠密的点云,可以加快后续点云处理的效率。
(3)点云统计滤波器
点云统计滤波器,是指根据统计学的方法,判断点云的密集程度,点云密集程度低的点一般是测量噪声导致的杂散点,然后根据点云的密集程度可以删除点云杂散点,在点云层面减少测量过程中噪声的影响。由于测量环境或多或少会有些噪声影响,统计滤波器对大多数的三维测量数据都有效。
(4)点云半径滤波器
点云半径滤波器,是指在点云中某个点的邻域半径区间内,计算区间内所有点的数量,当点的数量少于某个阈值时删除该点。半径滤波器相当于对点云数据进行图像腐蚀一样的操作,减少点云的数量,突出点云的密集性,处理速度较快,在某些点云处理过程中需要用到。
下面为对点云数据的具体处理过程。
首先,通过统计滤波器去除点云中的离散点。
其次,由于面阵三维测量的点云数据很密集,影响点云配准的速度,因此需要对点云数据进行降采样,使用体素格滤波器进行点云降采样。
接着,为了提高点云配准的准确度,在点云配准之前,对两组点云数据使用直通滤波器进行处理。
最后,两组点云数据一般是相邻的两组测量数据,并且机器人能够读取机器人末端位置信息,根据相邻两组测量数据的前后左右关系,通过直通滤波器截取保留两组点云的重合部分。对两组点云的重合部分进行点云配准,并且得到点云配准参数,再通过该点云配准参数对直通滤波器处理之前的两组完整点云数据进行配准。
本发明的点云配准使用ICP(Iterative Closest Point)算法,其又称为迭代最近点法,是指寻求最佳的点云P和点云Q之间的旋转平移变换关系,使得这两组点云能够合并成一组完整正确的点云数据。ICP算法的基本原理流程图如图5所示。
在图5的ICP算法原理流程图中,包括如下步骤:
步骤a:从点云数据P中采集出部分点集p i ,所述点集p i 的采集方法包括均匀采集、随机采集和法矢采集。
步骤b:采集出点云数据Q中的对应点集q i ,所述点集q i 的采集方法包括点到点、点到投影和点到面;
步骤c:根据点集p i 和点集q i 计算出点集变换矩阵M pq ,使得误差函数最小。其中,变换矩阵M pq 一般是通过四元数法和SVD奇异值分解法求解出,包含旋转矩阵R和平移矩阵t。误差函数E(R,t)表达式如下:
Figure 217633DEST_PATH_IMAGE008
(7)
步骤d:判断平均距离d是否小于预设阈值,
Figure 745566DEST_PATH_IMAGE009
,平均距离d的计算公式如 下:
Figure 4509DEST_PATH_IMAGE010
(8)
对于ICP算法,有较多优点,也有相应不足。ICP算法优点:可以获得非常精确的配准效果;不必对处理的点云进行分割和特征提取;在较好的初值情况下,可以得到很好的算法收敛性。ICP算法的不足之处:在搜索对应点的过程中,计算量非常大;在寻找对应点集时,可能会产生一定数量的错误对应点,影响最终的匹配结果。
本发明的基于点云滤波的ICP点云精配准方法,是针对面阵三维测量点云数据的特点,通过对点云数据进行预处理,以实现更加快速、稳定和高精度的点云精配准。本发明对加了点云滤波的ICP点云精配准和没加点云滤波的点云精配准方法做了相关对比实验,下面给出其具体对比分析。
首先点云滤波的有无,对ICP点云配准的时间上有很大差别。对五组点云数据进行了对比实验,每组点云数据包含A点云和B点云,即对A点云和B点云实现ICP点云精配准。不采用点云滤波的ICP点云精配准耗费时间情况如表1所示。
表1不采用点云滤波的ICP点云精配准耗时情况
耗时情况 第一组点云 第二组点云 第三组点云 第四组点云 第五组点云
A点云数量 340550 351820 407848 332053 791027
B点云数量 305404 449810 258202 418224 658613
耗费时间/秒 106.470 76.277 65.543 66.547 54.378
ICP点云配准的一个缺点就是非常耗时,本节使用基于点云滤波的ICP点云配准方法,能够有效减少ICP点云配准的时间。同样对上述五组点云数据进行点云精配准,采用了点云滤波的ICP点云配准耗费时间情况如表2所示。
表2采用了点云滤波的ICP点云精配准耗时情况
耗时情况 第一组点云 第二组点云 第三组点云 第四组点云 第五组点云
A点云数量 21222 23095 25477 30030 42216
B点云数量 21530 26824 19060 32112 37660
耗费时间/秒 1.802 0.603 0.863 1.128 0.873
其次对于点云差异较大且点云公共区域较小时,只使用ICP点云配准效果不稳定,容易出现误匹配。因此,使用基于点云滤波的ICP点云配准方法,通过截取A点云和B点云的公共区域进行ICP点云配准,然后再用该配准参数对A点云和B点云进行全局配准,能够实现精度较高且稳定的ICP点云配准。通过截取公共区域的匹配情况如下,图6(a)是配准前的点云图,点云差距较大且公共区域较小。图6(b)是没有使用点云滤波的ICP配准效果,匹配后还存在点云偏差。图6(c)是截取公共区域后的ICP配准情况,由于没有其余部分点云的干扰,点云配准情况较好。图6(d)是使用了图6(c)图中的匹配结果参数,对原始点云的匹配结果,匹配效果和图6(c)图中的匹配效果相同,效果较好。因此,基于点云滤波的ICP点云配准方法,不仅能够有效地提高点云配准的时间效率,同时也能提高点云配准的稳定性和精度。
有益效果:
(1)双目加投影仪的具有较高面阵三维重建精度;
(2)手眼标定粗拼接后存在微小偏差,点云配准精拼接后,能够矫正偏差,得到完整的物体三维形貌点云数据;
(3)使用双目相机矫正后的手眼标定矩阵能够减小粗拼接的点云偏差;
(4)机器人平移运动下,手眼标定粗拼接的效果较好,不使用点云配准精拼接方法也能得到较好的拼接效果;当机器人变姿态运动时,手眼标定粗拼接会存在微小偏差,此时需要ICP点云精配准的方法来矫正点云偏差。基于点云滤波的ICP点云配准方法,不仅能够有效地提高点云配准的时间效率,同时也能提高点云配准的稳定性和精度
综上,本发明的负责大构件机器人面扫描三维重建方法能够实现高速度和高精度的大型物体三维形貌测量。
附图说明
图1 投射的7位格雷码编码图案;
图2基于格雷码的相位展开图;
图3 中值滤波前后效果:(a)具有毛刺误差的相位展开图;(b)中值滤波器消除毛刺的相位展开图;
图4基于点云滤波的ICP点云配准方法流程图;
图5 ICP点云配准算法原理流程图;
图6点云配准对比效果:(a)配准前的点云图;(b)没有使用点云滤波的ICP配准效果;(c)截取公共区域后的ICP配准情况;(d)基于点云滤波的ICP点云配准效果;
图7 双目标定实验过程:(a)棋盘格标定靶图片;(b)其中一个相机的双目标定图片;
图8 单位面阵三维扫描数据的精度验证:(a)球面点云图;(b)球面拟合半径误差图。
图9总体拼接效果:(a)扫描物体实物;(b)物体左侧面阵三维点云;(c)物体右侧面阵三维点云;(d)拼接前的两组原始点云;(e)手眼标定粗拼接效果;(f)点云配准精拼接效果;
图10 双目相机协同手眼标定的效果验证:(a)相机1手眼标定矩阵粗拼接效果;(b)矫正后的手眼标定矩阵粗拼接效果;
图11机器人平移运动和变姿态运动的手眼标定粗拼接区别:(a)光头强玩具模型实物图;(b)平移运动手眼标定粗拼接效果图;(c)变姿态运动手眼标定粗拼接效果图;
图12大构件的面阵三维扫描效果:(a)大尺寸的键盘实物图;(b)两次面阵三维数据拼接效果图;(c)总体12次面阵三维数据拼接效果叠加图;(d)拼接完成后的键盘表面完整点云效果图。
具体实施方式
下面结合附图和具体实施例,进一步阐明本发明,应理解这些实施例仅用于说明本发明而不用于限制本发明的范围,在阅读了本发明之后,本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
实施例1:
基于面结构光与工业机器人的大构建三维形貌测量系统,由两部分构成,即工业机器人系统和单位面阵三维测量系统。工业机器人为安川ERER-MA02010-A00-C型工业机器人。单位面阵三维测量系统,由两台相机、镜头和一台投影仪组成。相机使用的是两个BasleracA1920-40gm型号工业相机,镜头均为12mm镜头。投影仪为TI公司DLP Lightcrafter 4500型投影仪,投影分辨率为912×1140,投影波段为420-720nm,工作范围为0.5-2m。相机和投影仪都支持C++编程控制,且都支持外部触发。单位面阵三维测量系统通过钳子等夹具,能够固定在工业机器人的末端,从而形成大构件三维形貌测量系统。
系统搭建之后,随后开始系统的标定。整个测量系统需要做两个标定,即双目标定和手眼标定。由于两个标定都使用棋盘格标定靶,因此实验中把双目标定和手眼标定的实验过程同步进行,做一组双目相机标定的实验,同时得到手眼标定需要的数据。本发明使用的棋盘格标定靶(如图7(a)所示)和相机拍摄的标定图片(如图7(b)所示)。
双目标定和手眼标定的实验过程如下:棋盘格放在机器人工作台上面全程固定不动;变换一次机器人位姿,待机器人静止不动后,双目相机分别对棋盘格标定靶成像并保存拍摄的图片,同时记录下机器人在该位姿下的末端坐标系到基底坐标系的转换关系参数,双目相机图片和机器人参数作为一组标定数据;记录完一组标定数据后,再次变换机器人位姿,待机器人停止后,重复上一步的操作,记录下另一组标定数据;然后重复上述实验步骤,直至获取二十多组的标定数据。
实验获取了标定数据之后,首先需要标定出双目相机的参数。对于相机标定,本节使用张正友相机标定方法。通过双目标定后,能够得到双目相机的标定参数,标定出来的参数如表3所示。
表3 双目相机的标定参数
Figure 696522DEST_PATH_IMAGE012
另外,需要进行手眼标定。双目标定能够得到标定出两个相机各自的标定参数,然后利用每个相机单独的相机参数标定出各自相机的手眼转换矩阵。具体的手眼标定步骤如下:
(1)分离双目标定参数,得到两个相机各自的相机参数,以及相机拍摄的20组外参矩阵;
(2)按照相机标定的棋盘格图片顺序,将记录的20组机器人数据与相机外参对应起来;
(3)根据前文介绍的基于棋盘格的手眼标定方法原理,通过相机外参矩阵和机器人参数计算出两个相机各自的手眼标定矩阵;
(4)然后通过双目转换矩阵和相机2的手眼转换矩阵,对相机1的手眼转换矩阵进行矫正;
(5)矫正后的相机1的手眼转换矩阵,即为测量系统中待求的手眼转换矩阵。
手眼标定过程中,直接求解出的相机1和相机2的手眼转换矩阵如下:
Figure 886063DEST_PATH_IMAGE014
(10)
Figure 347132DEST_PATH_IMAGE016
(11)
通过相机2的手眼标定矩阵X 2 和双目转换矩阵H 12 结合后,得到新的相机1手眼转换矩阵如下:
Figure 721785DEST_PATH_IMAGE018
(12)
然后,根据
Figure DEST_PATH_IMAGE019
矩阵对相机1的手眼标定矩阵进行矫正,得到最终的测量系统手眼标定转 换矩阵如下:
Figure 88044DEST_PATH_IMAGE021
(13)
经过测量系统的搭建和标定之后,就能使用本节的大构件三维形貌测量系统对实际物体进行三维扫描测量。整体系统标定完成之后,本节进行了以下五组相关验证实验。
(1)首先为验证单位面阵三维扫描数据的准确性,使用标准球进行格雷码结合相移的单位三维面阵测量精度验证实验。不移动机器人,只是通过双目相机和投影仪对标准球进行单次重建,然后对重建出来的标准球点云数据进行球面拟合。拟合的球面点云图(如图8(a)所示)和球面拟合半径误差(如图8(b)所示),球面上所有点到拟合球心点的半径距离均方根误差为0.082392mm,说明双目加投影仪的具有较高面阵三维重建精度。
(2)然后使用该测量系统从两个不同的侧面扫描了某娃娃玩具模型,用来验证手眼标定粗拼接和点云配准精拼接的效果,娃娃玩具模型的实物图如图9(a)所示。从两个侧面扫描物体,需要机器人带动相机和投影仪运动到不同的姿态和位置,从两个侧面扫描的三维点云数据分别如图9(b)和9(c)所示。扫描出两个单位面阵点云数据之后,拼接前的两组原始单位点云的叠加图如图9(d)所示。单位扫描完成后,使用手眼标定粗拼接方式进行三维拼接,手眼标定粗拼接的结果存在微小偏差,粗拼接效果图如图9(e)所示。手眼标定粗拼接完成后,使用基于点云滤波的ICP点云配准方法,对粗拼接结果进行精拼接,纠正粗拼接存在的微小偏差,点云配准精拼接效果图如图9(f)所示。由图9(f)的拼接效果可知,手眼标定粗拼接后存在微小偏差,点云配准精拼接后,能够矫正偏差,得到完整的物体三维形貌点云数据。
(3)其次对于双目相机协同进行手眼标定的方法,对单相机的手眼标定矩阵和双目相机矫正后的手眼标定矩阵拼接效果进行验证。由于获取的单位面阵三维测量数据是基于相机1的相机坐标系,只用相机1的手眼标定矩阵能够实现手眼标定粗拼接。同样使用上图的娃娃玩具模型,变姿态的扫描物体的左右两侧,展示点云侧面的拼接效果,相机1手眼标定矩阵粗拼接效果如图10(a)所示。另外,利用双目相机转换矩阵和相机2的手眼标定矩阵,能够对相机1的手眼标定矩阵进行矫正,双目相机矫正后的手眼标定矩阵粗拼接效果如图10(b)所示。由图10(b)的效果可知,只使用相机1手眼标定矩阵的粗拼接效果,点云偏差相对较大,使用双目相机矫正后的手眼标定矩阵能够减小粗拼接的点云偏差。
(4)另外对于机器臂平移运动和变姿态运动,验证大构件三维形貌扫描的差别。一般的工业机器人会有两种运动方式,即变姿态运动和平移运动。变姿态运动是指只改变机器人的六轴角度,使机器臂末端运动,这时机器臂末端是旋转运动且不规则。平移运动是指控制机器臂末端,按照机器人基底坐标系的三轴方向进行运动,此时机器臂末端平移运动且沿着坐标系的三轴方向。实验中发现,机器人单纯做平移运动时,只使用手眼标定粗拼接也能达到很好的效果;当机器人做变姿态运动时,只使用手眼标定粗拼接的效果会存在微小偏差,这时必须使用点云配准精拼接的方法来矫正点云偏差。为了验证机器人平移运动和变姿态运动的差别,本节扫描了某光头强玩具模型,该光头强玩具模型的实物如图11(a)所示。机器人做平移运动时,扫描了两次单位面阵三维测量数据,只使用手眼标定粗拼接的效果图如图11(b)所示,颜色不同的点云代表不同的单位面阵三维测量数据。机器人做变姿态运动后,也同样扫描了两次光头强玩具模型,使用手眼标定粗拼接的效果如图11(c)所示。由效果图可知,机器人平移运动下,手眼标定粗拼接的效果较好,不使用点云配准精拼接方法也能得到较好的拼接效果;当机器人变姿态运动时,手眼标定粗拼接会存在微小偏差,此时需要ICP点云精配准的方法来矫正点云偏差。
(5)最后使用测量系统扫描了某大型键盘,键盘尺寸是45×23cm,一般的面阵三维测量方法很难一次性扫描出该键盘的完整三维形貌数据,因此使用本节的大构件三维形貌测量方法来获取该键盘的表面完整三维点云数据。该键盘的实物图如图12(a)所示。由于本节使用的面阵三维测量方法重建的单位面积较小,为获取键盘表面完整的三维点云数据,实验过程中使用该测量系统对键盘进行12次单位面阵扫描,在整体拼接的过程中,相邻两组的单位面阵三维数据拼接效果如图12(b)所示。然后,对这12组单位面阵三维测量数据进行手眼标定粗拼接,再对相邻点云之间进行点云配准精拼接,拼接后的12组点云数据能够组成完整的键盘表面三维形貌,拼接后的12组点云图叠加效果如图12(c)所示,其中不同颜色的点云代表着不同的单位面阵三维测量数据。把拼接后的12组点云数据放在一起,就能够组成键盘表面完整的点云数据,完整的键盘三维形貌点云图如图12(d)所示。由实验拼接效果可知,本节的基于工业机器人的大构件三维形貌测量系统能够实现高精度的大型物体三维形貌测量。

Claims (6)

1.一种复杂大构件机器人面扫描三维重建方法,其特征在于:包括
步骤1:采用格雷码编码图像结合相移的结构光编码方法进行面阵三维形貌测量;
步骤2:将面阵三维形貌测量系统和工业机器人结合起来,进行系统标定;
步骤3:对采集的点云数据进行处理,进行粗拼接和精拼接处理;
所述步骤1中的面阵三维形貌测量包括依序的双面相机标定、光栅投影与图像同步采集、编码图像的相位展开、基于相位的双目立体匹配和空间点的三维坐标求解;
所述步骤2中的系统标定为双目标定和手眼标定结合;
所述步骤3所述粗拼接为使用手眼标定方法,把工业机器人工作区域内的三维测量数据统一工业机器人基底坐标系;所述精拼接为通过点云处理的方法,把粗拼接存在的微小误差消除掉。
2.根据权利要求1所述的复杂大构件机器人面扫描三维重建方法,其特征在于:所述双面相机标定包括根据相机模型标定出两个相机的内参和外参,确定相机各自的相机坐标系和世界坐标系的组成参数;标定出两个相机坐标系之间的转换关系,得到相应的参数,确定双目视觉测量模型;
所述光栅投影与图像同步采集包括将编码面结构光图案通过投影仪投射到物体表面;投影仪投影时,将正弦条纹图案和格雷码图案按照顺序依次投射出去;此时相机与投影仪同步,每当一幅编码图案投出去后,相机立刻拍摄被物体表面调制的编码图像;然后把图像保存下来,用于后续处理;
所述编码图像的相位展开包括对保存的图像进行解码,获取每幅图像像素点的相位值;
所述基于相位的双目立体匹配包括在编码图像经过解相位之后,左右图像得到每个像素点的相位值;同时根据双目立体视觉原理,遍历图像上的所有像素点,求出左右图像上的所有极线对;在极线对的两条极线上,搜索相位值相同的点,获得左右图像上的匹配点对;
所述空间点的三维坐标求解包括根据物点、像点、两个相机相对位置关系以及双目标定的结果参数,在图像二维坐标恢复物体三维坐标下,通过最小二乘法求解物点的三维坐标值。
3.根据权利要求2所述的复杂大构件机器人面扫描三维重建方法,其特征在于:所述编码图像的相位展开使用格雷码算法:通过将格雷码编码图像加载到投影仪内,经投影仪投射后被相机捕获,条纹解码的过程则通过将所采集的灰度图像经灰度映射外加设定阈值完成;如果某一点像素解码为B(x,y),则该点的绝对相位值为:
Figure 580907DEST_PATH_IMAGE001
(1)
假设包裹相位的分布函数为φ 1 (x,y),φ(x,y)表示经过物体表面调制后的相位分布,采用格雷码直接解码获得绝对相位的分布函数为Φ 1 (x,y),经中值滤波后所获取的绝对相位分布为Φ 2 (x,y),则理想绝对相位Φ(x,y)可通过如下公式获得:
Figure 837445DEST_PATH_IMAGE002
(2)
公式(2)中,round()为四舍五入函数,最终所获得理想相位便是解码准确时所获得绝对相位。
4.根据权利要求2所述的复杂大构件机器人面扫描三维重建方法,其特征在于:所述双目标定和手眼标定结合包括:
首先:硬件设定:设置棋盘格相对于机器人的底座固定不动,且始终位于双目相机的公共视场内,相机捕获棋盘格图像;
其次:双目标定:控制工业机器人带动固定在末端的双目相机,以不同的位置姿态拍摄棋盘格的标定图片;双目相机标定获得到各自相机的内参、外参、两个相机之间的转换矩阵和相机1到相机2的转换矩阵记为M
最后:手眼标定,在进行手眼标定时,把双目相机看成两个没有联系的相机,再分别对各自单独的相机进行手眼标定;棋盘格手眼标定方法是针对单目相机,得出两个单独相机的各自手眼标定参数,相机1的手眼转换矩阵记为X 1 ,相机2的手眼转换矩阵记为X 2 ;对于一个相机1的相机坐标系下的三维点P 1 ,其在相机2的相机坐标系下的三维点P 2 和在机器人末端坐标系下的三维点P e ,有如下两组转换关系:
Figure 238470DEST_PATH_IMAGE003
(3)
Figure 5700DEST_PATH_IMAGE004
(4)
公式(3)中,先把P 1 点转到相机2的坐标系下,然后再相机2的手眼标定矩阵转换到机器人末端坐标系下;在公式(4)中,直接把P 1 点通过相机1的手眼标定矩阵转换到机器人末端坐标系下;由上面两个等式可得:
Figure 425049DEST_PATH_IMAGE005
(5)
在理想情况下,相机1和相机2的手眼标定矩阵存在如上关系;但在实际情况下,由于存在机器人绝对定位误差等原因,使得上式不完全相等,即存在手眼标定误差;因此,可以通过上式的约束条件来对手眼标定矩阵进行误差校正;在双目面阵三维测量过程中,一般把三维测量数据所在的三维坐标系定义为相机1的相机坐标系;由于把相机1的手眼标定结果作为三维拼接的转换矩阵,用相机2的手眼标定结果对相机1的手眼标定结果进行误差矫正;通过以下计算方式,能够得到更加准确的相机1手眼标定矩阵,如公式(6)所示:
Figure 603221DEST_PATH_IMAGE006
(6)
其中,X 0 即为通过相机2矫正后的相机1手眼标定矩阵。
5.根据权利要求1所述的复杂大构件机器人面扫描三维重建方法,其特征在于:所述精拼接选用基于点云滤波的ICP点云配准方法,包括:
首先:通过统计滤波器去除点云中的离散点;
其次:使用体素格滤波器进行点云降采样;
接着:根据相邻两组测量数据的前后左右关系,通过直通滤波器截取保留两组点云的重合部分;
最后:对两组点云的重合部分进行点云配准,并且得到点云配准参数,再通过该点云配准参数对直通滤波器处理之前的两组完整点云数据进行配准。
6.根据权利要求5所述的复杂大构件机器人面扫描三维重建方法,其特征在于:所述
点云配准包括步骤:
步骤a:从点云数据P中采集出部分点集p i ,所述点集p i 的采集方法包括均匀采集、随机采集和法矢采集;
步骤b:采集出点云数据Q中的对应点集q i ,所述点集q i 的采集方法包括点到点、点到投影和点到面;
步骤c:根据点集p i 和点集q i 计算出点集变换矩阵M pq ,使得误差函数最小;其中,变换矩阵M pq 一般是通过四元数法和SVD奇异值分解法求解出,包含旋转矩阵R和平移矩阵t,n代表点集个数;误差函数E(R,t)表达式如下:
Figure 507854DEST_PATH_IMAGE007
(7)
步骤d:判断平均距离d是否小于预设阈值,平均距离d的计算公式如下:
Figure 62463DEST_PATH_IMAGE008
(8)
其中
Figure 601898DEST_PATH_IMAGE009
CN202011114453.9A 2020-10-19 2020-10-19 一种复杂大构件机器人面扫描三维重建方法 Active CN112013792B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011114453.9A CN112013792B (zh) 2020-10-19 2020-10-19 一种复杂大构件机器人面扫描三维重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011114453.9A CN112013792B (zh) 2020-10-19 2020-10-19 一种复杂大构件机器人面扫描三维重建方法

Publications (2)

Publication Number Publication Date
CN112013792A true CN112013792A (zh) 2020-12-01
CN112013792B CN112013792B (zh) 2021-02-02

Family

ID=73527367

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011114453.9A Active CN112013792B (zh) 2020-10-19 2020-10-19 一种复杂大构件机器人面扫描三维重建方法

Country Status (1)

Country Link
CN (1) CN112013792B (zh)

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112710235A (zh) * 2020-12-21 2021-04-27 北京百度网讯科技有限公司 结构光测量传感器的标定方法和装置
CN112950765A (zh) * 2021-03-08 2021-06-11 北京理工大学 一种基于点云建模的腔体直线度检测方法
CN113012277A (zh) * 2021-02-03 2021-06-22 中国地质大学(武汉) 一种基于dlp面结构光多相机重建方法
CN113251926A (zh) * 2021-06-04 2021-08-13 山东捷瑞数字科技股份有限公司 一种不规则物体的尺寸测量方法及测量装置
CN113269674A (zh) * 2021-05-12 2021-08-17 深圳群宾精密工业有限公司 一种利用机器人坐标对固定3d图像空间自动拼接的方法
CN113421291A (zh) * 2021-07-16 2021-09-21 北京华睿盛德科技有限公司 利用点云配准技术和三维重建技术的工件位置找正方法
CN113505626A (zh) * 2021-03-15 2021-10-15 南京理工大学 一种快速三维指纹采集方法与系统
CN113532325A (zh) * 2021-06-08 2021-10-22 深圳市格灵精睿视觉有限公司 动态步数解相方法、电子设备及计算机可读存储介质
CN113643372A (zh) * 2021-10-18 2021-11-12 中国科学院自动化研究所 一种三维焊缝提取方法及系统
CN114862966A (zh) * 2022-05-10 2022-08-05 东北大学 一种基于一维码的多线阵相机快速标定带及标定系统
CN115790449A (zh) * 2023-01-06 2023-03-14 威海晶合数字矿山技术有限公司 一种狭长空间的三维形貌测量方法
WO2023060927A1 (zh) * 2021-10-14 2023-04-20 五邑大学 一种3d光栅检测方法、装置、计算机设备及可读存储介质
CN117333649A (zh) * 2023-10-25 2024-01-02 天津大学 一种动态扰动下高频线扫描稠密点云的优化方法
CN117333649B (zh) * 2023-10-25 2024-06-04 天津大学 一种动态扰动下高频线扫描稠密点云的优化方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103868460B (zh) * 2014-03-13 2016-10-05 桂林电子科技大学 基于视差优化算法的双目立体视觉自动测量方法
CN104408762A (zh) * 2014-10-30 2015-03-11 福州大学 利用单目和二维平台获取物体图像信息及三维模型的方法
CN109272523B (zh) * 2018-08-13 2021-03-16 西安交通大学 基于改进cvfh和crh特征的随机堆放活塞位姿估计方法
CN109345620B (zh) * 2018-08-13 2022-06-24 浙江大学 融合快速点特征直方图的改进icp待测物体点云拼接方法
CN110695982A (zh) * 2019-10-17 2020-01-17 南京隆越自动化科技有限公司 一种基于三维视觉的机械臂手眼标定方法和装置
CN111156925B (zh) * 2019-12-19 2021-12-28 南京理工大学 基于线结构光和工业机器人的大构件三维测量方法

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112710235A (zh) * 2020-12-21 2021-04-27 北京百度网讯科技有限公司 结构光测量传感器的标定方法和装置
CN113012277A (zh) * 2021-02-03 2021-06-22 中国地质大学(武汉) 一种基于dlp面结构光多相机重建方法
CN113012277B (zh) * 2021-02-03 2022-06-21 中国地质大学(武汉) 一种基于dlp面结构光多相机重建方法
CN112950765A (zh) * 2021-03-08 2021-06-11 北京理工大学 一种基于点云建模的腔体直线度检测方法
CN113505626A (zh) * 2021-03-15 2021-10-15 南京理工大学 一种快速三维指纹采集方法与系统
CN113269674A (zh) * 2021-05-12 2021-08-17 深圳群宾精密工业有限公司 一种利用机器人坐标对固定3d图像空间自动拼接的方法
CN113251926A (zh) * 2021-06-04 2021-08-13 山东捷瑞数字科技股份有限公司 一种不规则物体的尺寸测量方法及测量装置
CN113532325A (zh) * 2021-06-08 2021-10-22 深圳市格灵精睿视觉有限公司 动态步数解相方法、电子设备及计算机可读存储介质
CN113421291A (zh) * 2021-07-16 2021-09-21 北京华睿盛德科技有限公司 利用点云配准技术和三维重建技术的工件位置找正方法
CN113421291B (zh) * 2021-07-16 2023-10-24 北京华睿盛德科技有限公司 利用点云配准技术和三维重建技术的工件位置找正方法
WO2023060927A1 (zh) * 2021-10-14 2023-04-20 五邑大学 一种3d光栅检测方法、装置、计算机设备及可读存储介质
CN113643372A (zh) * 2021-10-18 2021-11-12 中国科学院自动化研究所 一种三维焊缝提取方法及系统
CN113643372B (zh) * 2021-10-18 2022-03-04 中国科学院自动化研究所 一种三维焊缝提取方法及系统
CN114862966A (zh) * 2022-05-10 2022-08-05 东北大学 一种基于一维码的多线阵相机快速标定带及标定系统
CN115790449A (zh) * 2023-01-06 2023-03-14 威海晶合数字矿山技术有限公司 一种狭长空间的三维形貌测量方法
CN117333649A (zh) * 2023-10-25 2024-01-02 天津大学 一种动态扰动下高频线扫描稠密点云的优化方法
CN117333649B (zh) * 2023-10-25 2024-06-04 天津大学 一种动态扰动下高频线扫描稠密点云的优化方法

Also Published As

Publication number Publication date
CN112013792B (zh) 2021-02-02

Similar Documents

Publication Publication Date Title
CN112013792B (zh) 一种复杂大构件机器人面扫描三维重建方法
CN110514143B (zh) 一种基于反射镜的条纹投影系统标定方法
CN109272570B (zh) 一种基于立体视觉数学模型的空间点三维坐标求解方法
CN112465912B (zh) 一种立体相机标定方法及装置
CN113160339B (zh) 一种基于沙姆定律的投影仪标定方法
Douxchamps et al. High-accuracy and robust localization of large control markers for geometric camera calibration
CN111981982B (zh) 一种基于加权sfm算法的多向合作靶标光学测量方法
WO2018201677A1 (zh) 基于光束平差的远心镜头三维成像系统的标定方法及装置
CN113920205B (zh) 一种非同轴相机的标定方法
Huang et al. A single-shot-per-pose camera-projector calibration system for imperfect planar targets
CN113841384B (zh) 校准装置,用于校准的图表和校准方法
CN111899290B (zh) 一种联合偏振和双目视觉的三维重建方法
CN115861445B (zh) 一种基于标定板三维点云的手眼标定方法
CN110599578A (zh) 一种真实感三维彩色纹理重建方法
CN114004901A (zh) 多相机标定方法、装置、终端设备及可读存储介质
CN114549660B (zh) 基于圆柱形自识别标记物的多相机标定方法、装置及设备
CN117450955A (zh) 基于空间环形特征的薄型物体三维测量方法
CN114993207B (zh) 基于双目测量系统的三维重建方法
CN113865514B (zh) 一种线结构光三维测量系统标定方法
CN114820798A (zh) 一种标定器匹配方法及装置
CN113066131A (zh) 一种相机内参的快速标定方法
Resch et al. Semi-automatic calibration of a projector-camera system using arbitrary objects with known geometry
Vera et al. Relevance of Accurately Determining the Center of Radial Distortion in 3-D Reconstruction
KR102584209B1 (ko) 오목 렌즈 배열을 이용한 집적 영상의 3차원 재구성 방법
CN117804381B (zh) 一种基于相机阵列聚焦结构光对物体三维重建方法

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220105

Address after: 210000 room 502-3, block a, Kechuang headquarters building, No. 320, pubin Road, Nanjing area, Nanjing Free Trade Zone, Jiangsu Province

Patentee after: Nanjing shangyuhe Technology Co.,Ltd.

Address before: 211505 building 527, Fukang Park, 399 Xiongzhou South Road, Longchi street, Liuhe District, Nanjing City, Jiangsu Province

Patentee before: Nanjing Zhipu Photoelectric Technology Co.,Ltd.

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220913

Address after: Room 409, building B, Xingzhi science and Technology Park, No. 6, Xingzhi Road, Nanjing Economic and Technological Development Zone, Nanjing, Jiangsu 210000

Patentee after: Nanjing nanxuan Heya Technology Co.,Ltd.

Address before: 210000 room 502-3, block a, Kechuang headquarters building, No. 320, pubin Road, Nanjing area, Nanjing Free Trade Zone, Jiangsu Province

Patentee before: Nanjing shangyuhe Technology Co.,Ltd.