Viewpoint Planning for Robot Photogrammetry Based on Initial Pose Estimation via Deep Learning
-
摘要: 针对机器人摄影测量中离线规划受初始位姿标定影响的问题, 提出融合初始位姿估计的机器人摄影测量系统视点规划方法. 首先构建基于YOLO (You only look once) 的深度学习网络估计被测对象3D包围盒, 利用PNP (Perspective-N-point)算法快速求解对象姿态; 然后随机生成机器人无奇异无碰撞的视点, 基于相机成像的2D-3D正逆性映射, 根据深度原则计算每个视角下目标可见性矩阵; 最后, 引入熵权法, 以最小化重建信息熵为目标建立优化模型, 并基于旅行商问题(Travelling saleman problem, TSP)模型规划机器人路径. 结果表明, 利用深度学习估计的平移误差低于5 mm, 角度误差低于2°. 考虑熵权的视点规划方法提高了摄影测量质量, 融合深度学习初始姿态的摄影测量系统提高了重建效率. 利用本算法对典型零件进行摄影测量质量和效率的验证, 均获得优异的位姿估计和重建效果. 提出的算法适用于实际工程应用, 尤其是快速稀疏摄影重建, 促进了工业摄影测量速度与自动化程度提升.Abstract: Aiming at the problem that offline planning of robot photogrammetry is affected by the initial pose calibration, a viewpoint planning method of robot photogrammetry system incorporating initial pose estimation is proposed. First, we construct a YOLO (you only look once)-based deep learning network to estimate the 3D bounding box of the measured object, and utilize the PNP (perspective-N-point) algorithm to quickly solve the object pose; Second, we randomly generate non-singular and collision-free viewpoints. Based on the 2D-3D forward and inverse mapping of camera imaging, we calculate the target visibility matrix under each perspective according to the depth principle; Finally, the entropy-weighted method is introduced, the optimization model is established with the goal of minimizing the reconstruction information entropy afterward the robot path is planned based on the TSP (travelling salesman problem) model. The results show that the translation error estimated via deep learning is less than 5 mm, and the angular error is less than 2°. The viewpoint planning method considering entropy weight improves the quality of photogrammetry. Simultaneously, the reconstruction speed is increased. It obtains excellent pose estimation and reconstruction results when utilizing the algorithm to verify the photogrammetric quality and efficiency of more typical parts. The proposed algorithm is extendable to practical engineering applications, especially for rapid sparse photogrammetry, improving the speed and automation of industrial photogrammetry.
-
Key words:
- Photogrammetry /
- robot /
- deep learning /
- viewpoint planning /
- visibility matrix /
- entropy weight method
-
传统的摄影测量技术是利用单相机从多个角度拍摄大量图片进行物体三维重建的技术[1], 在机器人系统[2-3]、遥感测绘、制造测量[4]、航空航天[5-6]等领域应用广泛. 近年来, 利用机器人进行测量检测的研究和应用逐渐增加, 视觉与机器人系统的集成提升了视觉测量与检测的效率和自动化水平. 然而, 机器人摄影测量系统中, 机器人相对测量对象初始位姿需要标定和人工输入, 导致机器人离线编程受到实际测量环境影响.
机器人摄影测量系统的核心问题是安装视觉系统后, 具有可视能力的机器人视点规划与路径规划问题. Kinnell等[7]对机器人三维视觉系统框架做出了详细论述. Kwon等[8]研究了机器人扫描系统, 提出重扫描策略, 缩短了扫描时间, 增强了机器人检测的稳定性. Raffaeli等[9]提出基于理论数模的机器人扫描路径离线编程算法, 但是离线程序依赖于在线校准. Li等[10]提出基于邻域搜索方法来确定机器人二维扫描的方向, 但其路径规划不连续, 需要后期多组测量数据的再次融合. 为实现以最小的拍摄视点获得最丰富的重建信息, 需要依据视点可见性对拍摄视点进行优化[11]. “Generate-test”方法在基于模型的视点规划中得到了广泛的应用[12]. 它生成大量的候选视点, 并计算每个视点下的可见矩阵. 然后将覆盖问题转化为可见矩阵的最优选择问题[12]. Jing等[13]采用迭代随机抽样和概率势场方法以非确定性方式生成候选视点. 在较高的覆盖率下, 保证了较好的视点规划效果. 然后, Jing等[14]将视点规划和路径规划结合为一个单一的规划问题, 取得较好的全局最优效果. 在实际测量中, 机器人视觉系统必须确定被测对象相对于自身的实际位姿, 并以此为初始姿态进行扫描路径的规划. 基于视觉的目标位姿测量方法主要有PNP (Perspective-N-point)[15]方法和三维刚体变换法[16], 这些方法需要精确已知被测对象特征, 提取并重建出必要的特征要素进行位姿计算, 虽然位姿计算精度高, 但是在可变化的环境中特征提取不稳定, 容易产生误差. 深度学习技术在位姿估计方面具有环境适应性强、自动化程度高等优势. 近年来, 利用卷积神经网络(Convolutional neural network, CNN)进行位姿估计的研究发展迅速[17-19]. PoseNet[20]网络直接利用CNN从单幅RGB图像中估计相机6D姿态, 但是其将旋转分量和平移分量单独输出, 训练过程必须调整参数来平衡两个分量的损失. 为避免这个问题, PoseCNN[21]网络则直接将位姿误差作为损失函数而更加适应于位姿优化. 文献[22-23]中, CNN没有直接用于位姿估计, 而是通过估计特征点然后计算三维姿态的方法, 提高了位姿估计的速度和精度. 二维物体识别与预测主要有Faster R-CNN (Faster region-CNN)、YOLO (You only look once)、SSD (Single shot multibox detector)[22]等网络结构, 且均有良好的预测效果. Tekin等[17]和Rad等[23]分别利用YOLO和SSD提出了基于单幅图像的位姿估计模型, 提升了位姿估计速度.
可见, 利用深度学习进行位姿估计, 并与机器人视点规划相结合, 是机器人自动化高效率摄影测量的新的方式. 本文针对机器人摄影测量系统中初始位姿估计、机器人扫描路径规划两个主要问题, 研究提升机器人摄影测量效率和自动化程度的新方法. 提出融合深度学习位姿估计的机器人摄影测量系统视点规划策略(第2节). 主要贡献在于: 1)基于YOLO模型与PNP算法实现了单相机快速位姿估计, 替代传统的初始位姿现场校准的工作(第2.1节); 2)提出基于相机成像2D-3D正逆性映射的被测模型可见性矩阵估计(第2.2节); 3)提出基于初始位姿约束和熵权法的机器人摄影测量视点规划和路径规划方法(第2.3节). 利用几种具有典型特征的零件进行实验验证, 获得了优异的位姿估计效果和摄影重建结果(第3节).
1. 测量系统组成及方法
摄影测量是利用多视角拍摄大量图像进行物体三维重建的技术. 机器人摄影测量系统组成包括: 机器人本体、视觉末端、计算机及控制器等. 视觉系统固定安装在机器人末端执行器, 视觉坐标系与机器人末端坐标系精确标定后, 视觉数据可变换到机器人坐标系下. 机器人摄影测量系统中, 机器人作为执行机构, 搭载相机达到指定位置, 并以优化的视角获取被测对象数据. 图1为机器人摄影测量系统简图.
如图1所示, 被测对象放置在固定的工作平台上, 机器人携带视觉系统从不同角度拍摄图像, 然后发送至摄影测量系统进行重建. 由于摆放位置随意性, 传统机器人摄影测量必须借助辅助测量设备标定机器人基坐标系{B}和物体坐标系{O}的相对变换${\boldsymbol{T}}_O^B={\boldsymbol{T}}_M^B{\boldsymbol{T}}_O^M$. 然后利用这个固定关系和视觉系统可见信息规划机器人末端位姿${\boldsymbol{T}}_E^B$. 这一过程依赖初始标定结果, 使得离线编程结果受到实际操作限制. 本文从逆向出发, 利用视觉系统自身可定位能力实现被测对象相对于视觉系统的变换矩阵${\boldsymbol{T}}_C^O$估计, 然后根据手眼矩阵${\boldsymbol{T}}_C^E$计算出机器人姿态, 并根据可见性规划机器人完成扫描的最优位姿.
如图2所示, 整个视点规划策略可分为3部分: 初始位姿估计、可见性矩阵估计、视点与轨迹规划. 被测对象随意放置在测量平台, 机器人视觉系统仅通过初始图像即获取被测对象位姿, 而无需标定. 获得初始姿态后再根据可见性信息规划视点和机器人路径.
视点规划策略具体如下: 1) 基于深度学习的初始位姿估计: 训练YOLO网络[17]获得网络权重. 被测对象放置在测量系统前时采集初始图像, 利用已训练的网络估计初始位姿. 2) 可见性矩阵计算: 随机生成测量视点. 根据相机成像的2D-3D正逆性映射, 以三角网格模型为对象, 根据深度原则判断该视角下每个三角面片可见性, 生成可见性矩阵. 3) 视点规划建模与求解: 引入熵权法计算该视角下被测对象可重建的信息权值, 约束初始视角, 优化求解视点规划模型和旅行商模型获得最优视点和机器人扫描轨迹.
2. 视点规划策略
2.1 基于深度学习的单幅图像初始位姿估计
采用深度学习方法以实现在可变测量环境下被测目标的稳定位姿估计[17]. 采用的网络基于YOLO模型, 其通过识别对象包围盒及质心, 并利用PNP方法计算对象位姿. 该网络可以实现单幅图像的快速位姿估计. 图3为实现单幅图像位姿估计的网络结构和流程.
如图3所示, 本文采用文献[17]的网络结构, 本网络中定义可信度函数为
$$c(x)=\left\{ \begin{aligned} &{{{\rm{e}}^{\alpha \left( {1 - \frac{{{D_T}(x)}}{{{d_{th}}}}} \right)}},}{\;\;\;{D_T}(x) < {d_{th}}}\\ &0,\;\;\;\;\;\qquad\qquad{\text{否则}} \end{aligned} \right.$$ (1) 其中, ${D_T}(x)$表示图像平面两个像素点的欧氏距离; $\alpha $为指数函数的控制参数, 根据实际训练经验设置为2; ${d_{th}}$为像素差异阈值. 根据指定的测量对象, 在训练过程中, 定义融合点误差和致信误差的复合损失函数, 不考虑被测对象种类识别的损失, 则
$$L={\lambda _{pt}}{L_{pt}} + {\lambda _{conf}}{L_{conf}}$$ (2) 其中, ${L_{pt}}$和${L_{conf}}$分别表示图像点和可信度的损失函数, 利用均方根误差计算; ${\lambda _{pt}}$表示坐标点损失系数, 设置其值为1; ${\lambda _{conf}}$为可信度损失系数, 训练时对图像的每个元胞进行搜索, 若该元胞内有目标物体, ${\lambda _{conf}}=5$, 否则${\lambda _{conf}}=0.1$.
本文利用ArUco[22]编码标志点作为全局控制点, 并以其中一个编码标志点作为坐标系原点建立坐标系. ArUco具有较稳定的可识别性, 且可实现自身的位姿计算. 真值位姿为相机坐标系相对于全局坐标系的变换矩阵. 利用真值位姿可计算被测对象相关的任意投影图像点. 即对于k个拍摄视角, 若测得编码坐标系相对相机坐标系的真值位姿为${{\boldsymbol{T}}_k}$, 目标物体三维包围盒角点及质心在该坐标系下的坐标为${{\boldsymbol{M}}_i}={({X_i},{Y_i},{Z_i},1)^{\rm{T}}},i=0,1,\cdots,8$, 相机参数矩阵为${\boldsymbol{K}}$, 则可根据相机成像原理获得图像点mi, 则
$$s{{\boldsymbol{m}}_i}={\boldsymbol{K}}{{\boldsymbol{T}}_k}{{\boldsymbol{M}}_i}$$ (3) 其中s为尺度因子. 获得图像点后, 根据图像分辨率大小对图像点进行归一化处理. 利用这些图像点拟合最小的2D包围盒, 并对2D包围盒的长宽进行归一化处理. 训练网络需要准备的真值数据包括: 被测目标3D包围盒8个角点的归一化坐标、质心归一化坐标、图像2D包围盒长宽归一化数值以及目标物体编码, 共21个参数. 根据被测对象不同, 预先准备训练数据包括: 原始图像(三通道彩色)、掩膜图像(二值图)、目标物体三维模型、相机参数以及上述21个数值. 测试网络的输入为原始图像, 其通过在随机生成的视点处拍摄得到, 输出为上述21个参数. 利用图3网络训练后进行目标特征点预测, 获得目标的包围盒图像坐标后, 利用PNP估计位姿. 记预测的被测目标包围盒角点及质心图像坐标为${{\tilde{\boldsymbol{m}}}_i}$, 其相机坐标系下归一化平面坐标${{\hat{\boldsymbol{m}}}_i}({u_i},{v_i},1)={{\boldsymbol{K}}^{ - 1}}{{\tilde{\boldsymbol{m}}}_i}$, 预测位姿${{\tilde{\boldsymbol{T}}}_k}({{\boldsymbol{R}}_k},{{\boldsymbol{t}}_k})$为3 × 4的矩阵, 根据相机3D-2D投影关系(式(3))和PNP算法, 建立位姿计算的方程式, 即
$$s\left[ {\begin{array}{*{20}{c}} {{u_i}} \\ {{v_i}} \\ 1 \end{array}} \right]=\left[ {\begin{array}{*{20}{c}} {{r_{11}}}&{{r_{12}}}&{{r_{13}}}&{{t_1}} \\ {{r_{21}}}&{{r_{22}}}&{{r_{23}}}&{{t_2}} \\ {{r_{31}}}&{{r_{32}}}&{{r_{33}}}&{{t_3}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{X_i}} \\ {{Y_i}} \\ {{Z_i}} \\ 1 \end{array}} \right]$$ (4) 其中, ${{\boldsymbol{R}}_k}=\left[ \begin{aligned} {{r_{11}}}\;\;{{r_{12}}}\;\;{{r_{13}}} \\ {{r_{21}}}\;\;{{r_{22}}}\;\;{{r_{23}}} \\ {{r_{31}}}\;\;{{r_{32}}}\;\;{{r_{33}}} \end{aligned} \right]$, ${{\boldsymbol{t}}_k}=\left[ \begin{aligned} {{t_1}} \\ {{t_2}} \\ {{t_3}} \end{aligned} \right]$. 消去s,
${u_i}=\dfrac{{{r_{11}}{X_i} + {r_{12}}{Y_i} + {r_{13}}{Z_i} + {t_1}}}{{{r_{31}}{X_i} + {r_{32}}{Y_i} + {r_{33}}{Z_i} + {t_3}}}$
${v_i}=\dfrac{{{r_{21}}{X_i} + {r_{22}}{Y_i} + {r_{23}}{Z_i} + {t_2}}}{{{r_{31}}{X_i} + {r_{32}}{Y_i} + {r_{33}}{Z_i} + {t_3}}}$
记${{\tilde{\boldsymbol{T}}}_k}={\left[ {{{\boldsymbol{r}}_1},{{\boldsymbol{r}}_2},{{\boldsymbol{r}}_3}} \right]^{\rm{T}}}$, 其中${{\boldsymbol{r}}_j}$为${{\tilde{\boldsymbol{T}}}_k}$第j (j = 1, 2, 3)行的转置向量. 则式(4)可表示为
$$\left\{ \begin{aligned} &{\boldsymbol{r}}_1^{\rm{T}}{{\boldsymbol{M}}_i} - {\boldsymbol{r}}_3^{\rm{T}}{{\boldsymbol{M}}_i}{u_i}=0 \\ & {\boldsymbol{r}}_2^{\rm{T}}{{\boldsymbol{M}}_i} - {\boldsymbol{r}}_3^{\rm{T}}{{\boldsymbol{M}}_i}{v_i}=0 \end{aligned} \right.$$ (5) 利用网络预测的8个角点坐标及质心坐标, 建立线性方程组
$$\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{M}}_0^{\rm{T}}}&0&{ - {u_0}{\boldsymbol{M}}_0^{\rm{T}}} \\ 0&{{\boldsymbol{M}}_0^{\rm{T}}}&{ - {v_0}{\boldsymbol{M}}_0^{\rm{T}}} \\ \vdots & \ddots & \vdots \\ {{\boldsymbol{M}}_8^{\rm{T}}}&0&{ - {u_8}{\boldsymbol{M}}_8^{\rm{T}}} \\ 0&{{\boldsymbol{M}}_8^{\rm{T}}}&{ - {v_8}{\boldsymbol{M}}_8^{\rm{T}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{r}}_1}} \\ {{{\boldsymbol{r}}_2}} \\ {{{\boldsymbol{r}}_3}} \end{array}} \right]={\bf{0}}$$ (6) 通过奇异值分解法对系数矩阵进行分解, 获得${{\tilde{\boldsymbol{T}}}_k}$的最小二乘解.
2.2 基于2D-3D映射的可见性矩阵计算
为优化视点, 以最少的图片数量完成对象重建, 建立以可见性矩阵为目标函数的优化模型. 可见性矩阵是描述在某个视角下对象表面某个部分是否可见的二值化矩阵. 本文研究三角网格模型, 可见性矩阵描述为
$${A_{ij}}=\left\{ \begin{aligned} &{1,\quad \text{三角面图}\;j\; \text{在视图} \;i\; \text{中可见}}\\ &{0,\quad \text{否则}} \end{aligned}\right.$$ (7) 本文提出一种基于相机模型的空间映射方法确定可见性矩阵.
首先, 在被测对象自身建立局部坐标系, 随机生成空间三维点, 观测方向法向于原点. 记初始视点为${\boldsymbol{V}}_1= (x_{1i},y_{1i},z_{1i},\delta_{1i},\varphi_{1i}, \gamma_{1i})$, 其中,
$$\left\{\begin{aligned} &{x_{1i}}={r_i}\cos {\alpha _i}\cos {\beta _i} \\ &{y_{1i}}={r_i}\cos {\alpha _i}\sin {\beta _i} \\ &{z_{1i}}={r_i}\sin {\alpha _i} \end{aligned}\right.$$ (8) $$\left\{ \begin{aligned} & {\delta _{1i}}=\arcsin \left(\sqrt {\frac{{x_{1i}^2 + z_{1i}^2}}{{x_{1i}^2 + y_{1i}^2 + z_{1i}^2}}}\right) \\ & {\varphi _{1i}}=\arctan \left(\frac{x_{1i}}{z_{1i}}\right) \\ &{\gamma _{1i}}=0 \end{aligned} \right.$$ (9) 其中, ${r_i},\; {\alpha_i},\; {\beta_i}$为随机生成, $\delta_{1i},\varphi_{1i}, \gamma_{1i} $为欧拉角. 按照式(8)和式(9)生成的随机视点分布在空间一个局部环体中, 与实际拍摄视点相符. 根据手眼矩阵将V1变换到机器人末端坐标系, 然后根据机器人逆运动学求解得到机器人关节角度. 若求解失败或关节角度超出机器人理论关节角范围, 则将该视点移除. 初步滤除后的视点记为${{\boldsymbol{V}}_2}$. 控制机器人到达${{\boldsymbol{V}}_2}$后拍摄图像, 利用第2.1节深度学习的方法从图像中提取被测物体并估计当前姿态, 记为${{\boldsymbol{V}}'_2}$. ${{\boldsymbol{V}}_2}$为设计的视点, 位姿信息为已知. ${{\boldsymbol{V}}'_2}$为估计的参数, 位姿信息是估计得到的.
视点${{\boldsymbol{V}}'_2}$中每一个视点对应一个变换矩阵${\boldsymbol{T}}'$, 其是从被测对象坐标系${{\boldsymbol{O}}_m}$到相机坐标系OC的变换. 若空间三角面片3个顶点在${{\boldsymbol{O}}_m}$下坐标为${{\boldsymbol{P}}_s}\;(s=1, 2, 3)$, 利用式(3)将其投影到图像平面得到点${{\boldsymbol{p}}_s}\;(s=1, 2, 3)$. 记空白图像为I0, 分辨率与相机分辨率相同, 像素值均为0. 在图像平面I0顺序连接${{\boldsymbol{p}}_s}$3个像素点, 填充3个点所包围的三角形区域, 得到图像I1. 查找I1中非零像素点集${{\boldsymbol{m}}_{no}}$. 对于任意点$({u_i},{v_i}) \in {{\boldsymbol{m}}_{no}}$, 根据相机模型的2D-3D映射原理, 可重构出其在相机坐标系下的空间三维点${\boldsymbol{M}}_{no}^i\;(X_{no}^i,Y_{no}^i,Z_{no}^i)$, 即求
$$C_{no}^i{\boldsymbol{M}}_{no}^i={\bf{0}}$$ (10) 其中, $C_{no}^i=\left[\begin{aligned} &{C_1}{{{\boldsymbol{T}}'}_{}} \\& {C_2} \end{aligned} \right]$, ${C_1}=\left[ \begin{aligned} { - {f_x}}\;\;\;0\;\;\;{{u_i} - {u_0}} \\ \;\;0\;\;\;{ - {f_y}}\;\;{{v_i} - {v_0}} \end{aligned} \right]$, C2 = [a, b, c, d], $({u_0},{v_0})$为图像主点, a, b, c, d为三角面片形成的空间平面的方程式系数. 根据式(10)便可求得图像I2中任意图像点对应的三角面片上对应空间点的三维坐标. 三角面片中任意空间点存在与其他三角面片深度交叉, 根据最近原则确定该编号的三角面片是否可见. 可见性矩阵求解算法如下:
算法 1. 可见性矩阵计算算法
输入. 候选视点集合${{\boldsymbol{V}}'_2}$; 模型的所有三角面片集合${\boldsymbol{P}}$; 相机参数矩阵${\boldsymbol{K}}$.
输出. 可见性矩阵${{\boldsymbol{A}}_{ij}}\;(i \in length ({{\boldsymbol{V}}_2)},$其中$length(\cdot) $表示集合的大小, $j \in {\boldsymbol{P}}) $.
1: for each ${v_i} \in {{\boldsymbol{V}}'_2}$ do
2: initializing image Ii = all(255), Ic = all(0), Aij = 0. struct Pg and downpixel value.
3: for each ${{\boldsymbol{P}}_j} \in {\boldsymbol{P}}$ do
4: initializing image Ij = all(0), Id = all(1.0).
5: project ${\boldsymbol{P}}_j^s$ onto image plane I0 with (3), get ${\boldsymbol{p}}_j^s$;
6: fill I0, Ii, with ${\boldsymbol{p}}_j^s$and downpixel, count none zero pixel index in Ii, get Nidx.
7: for each $n \in {N_{idx}}$, $({u_i},{v_i}) \in {{\boldsymbol{m}}_{no}}$ do
8: calculate depth d of (ui, vi) with (10);
9: if $d < {I_d}({u_i},{v_i})$ do
10: ${I_d}({u_i},{v_i}) $ = d;
${I_c}({u_i},{v_i}) $ = $255-downpixel $
11: end if
12: end for
13: ${I_i}={I_i} - {I_j}$;
14: add j and Nidx into Pg
15: end for
16: for each $k \in $ Pg do
17: if size of $N_{idx}^k < 5$ continue;
18: if
${I_i}(mid(N_{idx}^k))$ $<255-2\times downpixel $
continue;
19: let $A(i,N_{idx}^k.idx)=1$
20: end for
21: end for
22: return Aij
2.3 测量视点优化与轨迹规划
为实现以最小图片数量对被测对象进行完整重建, 本文以最小化视点数量为目标函数, 可见性为约束条件, 建立基于最小化信息熵的优化模型$, $ 即
$$ \begin{split} &f=\min \sum\limits_{i \in {{{\boldsymbol{V}}}'_2}} {{w_i}{x_i}} \\ &{\rm{s.t}}. \left\{ \begin{aligned} &{x_1}=1 \\ &Ax \geq \gamma \\ &{x_i}({x_i} - 1)=0 \end{aligned} \right. \end{split} $$ (11) 不失一般性, ${x_1}=1$表示第一个视点为必须的视点, 其是当被测对象放置于测量系统前第一次拍摄的视点. 同时, 第一个视角下的位姿是利用深度学习方法估计得到的; $\gamma$是一个向量, 表示每个三角面片需要被重复拍摄的次数, 在摄影测量中, 根据重复测量的特征点计算相对位姿, 故必须保证每个特征点至少被测量2次, 本文设定为4; ${w_i}$表示视点${x_i}$下被测物体能够被重建的权重, 本文利用信息熵权的方法确定. 根据第2.3节算法流程可得每个视角下无遮挡的投影图像${I_c}$, 其对应一个投影矩阵${\boldsymbol{T}}({I_c})=\left[ {{\boldsymbol{R}}({I_c})|{\boldsymbol{t}}({I_c})} \right]$, 考虑到拍摄视角距离被测物体越近意味着更有可能实现特征识别和重建, 所以建立新的投影矩阵${\boldsymbol{T}}'({I'_c})=\left[ {{\boldsymbol{R}}({{I'}_c})|{\boldsymbol{t}}'({{I}'_c})} \right]$, 其中, ${\boldsymbol{R}}({I'_c})={\boldsymbol{R}}({I_c})$, ${\boldsymbol{t}}'({I'_c})$表示相机的最近可测深度, 其与相机景深有关. 利用${\boldsymbol{T}}'$将被测物体投影可得图像${I'_c}$. 定义${v_i}$个视点可重建概率
$$P({v_i})=\frac{{size({I_c})}}{{size({{I}'_c})}}$$ (12) 其中, size(I)表示图像I中像素值非零元素的个数. 概率越大, 该视点被选作最终视点的可能性越大. 本文目标函数为最小化, 根据信息熵理论, 定义信息熵为
$${w_i}= - \lg (P({v_i}))$$ (13) 信息熵与概率成反函数, 所以信息熵越小, 目标函数越小. 信息熵表明该视角下能够重建的信息值, 可重建信息熵值越小, 获得的优化的视点数越小. 利用线性规划的方法求解式(11), 可得到优化的视点集合V3.
下一步即是根据最短路径原则, 为机器人规划运动轨迹. 从${v_i}$到${v_j}$有多条路径, 定义机器人在两个视点之间的路径长度$d({v_i},{v_j})$为各个结点间距离总和. 为避免机器人运动至奇异点, 通常采用插值的方法在两个视点之间插入中间点, 本文采取自适应随机游走算法(Adaptive random walks, ARW)[24]对插值点和整个运动轨迹进行优化计算, 利用Robwork机器人编程库实现[25]. 可建立机器人路径规划的旅行商问题(Travelling saleman problem, TSP)模型
$$ \begin{split} &f=\min d({v_i},{v_j}){x_{ij}}\\ &{\rm{s.t}}.\left\{ \begin{aligned} &\sum\limits_{{v_i} \in V,{v_i} \ne {v_j}} {{x_{ij}}} =1 \\ &\sum\limits_{{v_j} \in V,{v_i} \ne {v_j}} {{x_{ij}}} =1 \\ &{x_{ij}}({x_{ij}} - 1)=0 \end{aligned}\right. \end{split} $$ (14) 利用MATLAB中求解TSP的工具箱进行求解, 可以得到最佳扫描路径.
3. 实验与结果分析
3.1 单幅图像位姿估计结果分析
本文选择一个3D打印件作为算法验证对象. 在实验中, 根据式(8)和式(9)生成1600个初始视点, 利用机器人逆向运动学和关节角度范围进行初步滤除, 获得1443个有效视点. 从中随机选择1155个视点拍摄图像进行训练, 其余288个视点获得的图像进行测试, 图像分辨率为1024$\times $768像素. 根据第2.1节准备测试数据进行网络训练. 目标物体与空间标志点的相对位置通过文献[15]方法求解. 训练损失与精度变化曲线如图4(a)所示, 测试集平移和旋转误差如图4(b)所示.
其中测试精度是每次训练后, 预测的包围盒角点与真实角点像素偏差小于5 pixels的角点所占的比率. 由图4(a)可知, 测试精度最终稳定在100%, 表明所有预测的目标包围盒均满足像素偏差小于5像素的要求. 平移精度和旋转精度定义如下:
$$\left\{\begin{aligned} &{E_{Trs}}=\frac{1}{N}\sum\limits_{i=1}^N {{{\left\| {Tr{s_{ie}} - Tr{s_{i0}}} \right\|}^2}} \\ &{E_{Ang}}=\frac{1}{N}\sum\limits_{i=1}^N {{{\left\| {An{g_{ie}} - An{g_{i0}}} \right\|}^2}} \end{aligned}\right. $$ (15) 其中, $Tr{s_{ie}}$, $Tr{s_{i0}}$分别表示估计的平移向量和真实平移向量; $An{g_{ie}}$, $An{g_{i0}}$分别表示估计欧拉角度和真实的欧拉角度. 从图4(b)可以看出, 训练稳定时, 平移精度为0.00481 m, 角度精度为1.6478°. 位姿估计的可视化结果如图5所示.
图5表示利用真值位姿和估计位姿将空间点重投影到图像平面的结果, 从图像可以看出, 重投影的空间物体与实际物体吻合得很好, 结合图4可得, 在机器人摄影系统中, 本网络结构可以实现较高精度的位姿估计. 单幅图像测得的空间位姿, 为机器人扫描路径规划提供准确的初始姿态.
3.2 可见性矩阵实验结果分析
本文以三角网格为研究对象, 首先将网格划分更加密集, 以提高可见性矩阵有效性和解析能力. 将测试集的288个视点作为候选视点, 按照第2.2节可见性矩阵的计算算法, 求解288个视点下的可见性. 部分计算结果如图6所示.
图6(a)表示直接将stl模型投影到图像平面所形成的投影图. 若图像中某个像素点只由一个空间点投影产生, 则该点的灰度值减少一次, 本文规定一次减少的灰度值为30. 则可知减少的次数越多, 灰度值越小, 表明该像素点由多个空间点投影产生, 即产生了将不可见空间面片投影到图像平面的错误投影. 由图6(a)可以看出, 投影区域的灰度值并不统一, 尤其在红色曲线内, 灰度值较小. 图6(b)是利用本文算法获得的投影图像, 灰度值完全相同, 表明所有像素点均是可见面片的投影. 图6(c)为对应视角下可见的三角面片, 其与投影图像相符. 图6(d)展示了该视角下拍摄图像, 然后进行位姿估计的结果, 这说明可以利用估计的位姿代替真值位姿进行可见性计算, 这更适用于实际工程应用. 通过标记每一个像素点对应空间面片的编号, 可以得到该面片在该视角下是否可见, 最终获得的可见性矩阵为二值化矩阵.
3.3 信息熵权法实验分析
获得可见性矩阵后, 利用式(11) ~ (14)即可以确定优化的测量视点和对应的扫描路径. 为探究本文提出的信息熵权法对于全局重建的影响, 本文对比考虑有权重和无权重情况下三维重建的质量, 结果如表1所示.
表 1 信息熵权有效性验证表Table 1 Effectiveness test for entropy weight对比实验 目标函数 优化的视
点数(个)点云数(个) A $x^*=\min \displaystyle \sum\limits_{i=1}^N { {w_i}{x_i} }$ 21 21509 B $x^*=\min \displaystyle \sum\limits_{i=1}^N { {x_i} }$ 20 18360 C $x^*=\min \displaystyle \sum\limits_{i=1}^N { {x_i} }$ 21 = 20 (B) +
1 (288)15344 表1中A、B、C三组实验, 分别设定了不同的目标函数. A组实验目标函数考虑了权重, 而B、C组实验均没有考虑权重. 在表1中, 优化的视点数表示经过优化求解后获得优化的视点数量, 点云数表示摄影测量系统重建的点云包含三维点的数量, 点云数越多表明被重建的特征越多, 重建质量越高. 对比A、B实验可看出, 增加了信息熵权的方法求解的视点数为21, 其视点数相对没有考虑权重的方法多了一个视点, 点云数增加了3149, 说明考虑权重后点云质量提高. C的目标函数和B相同, 所以得到的优化的视点数和B相同, 均为20个. 但是为说明输入图像数量不同对于重建质量的影响, 本文从候选的288个视点中随机选择一个视点, 故C共有21个视点. 比较重建点云数量可看出, 其重建的点云数量小于B, 说明同样21个视点情况下, 本文信息熵的方法可获得最优的视点数和三维重建质量.
进一步, 综合考虑初始姿态和权重情况下, 探究重建质量. 为提高重建速度, 在式(11)中控制每一个面片可见次数为2. 实验结果如表2所示.
表 2 综合权重和初始姿态下重建质量对比Table 2 Comparison of reconstruction quality with weight and first-sight pose有权重 无权重 有初始位姿约束 无初始位姿约束 有初始位姿约束 无初始位姿约束 视点索引 1, 13, 100, 113, 143,
149, 173, 189, 190,
196, 207, 269, 272, 28013, 100, 113, 143,
149, 173, 189, 190, 196,
207, 269, 272, 2801, 17, 28, 38, 45, 61,
66, 74, 89, 91, 92, 107,
113, 127, 185, 189,
207, 249, 269, 28014, 35, 43, 45, 56, 59,
73, 75, 89, 111, 127, 149,
162, 185, 189, 207, 249,
256, 274, 281三维点云 
点数量 10584 11451 7703 9571 注: 有初始位姿约束下，索引为1的视点需被约束保留. 表2对比了考虑权重和初始姿态情况下重建质量. 不失一般性, 统一将索引为1的视点作为初始视点. 其中有初始位姿约束表示约束候选视点中索引为1的视点为必选视点, 然后求解最优视点; 无初始位姿约束表示不考虑视点1是否为必选而直接求解最优视点. 表2结果保证了只有初始位姿是不同的. 对比表2第2列和第3列, 考虑初始姿态时, 视点数量增加1个, 但重建的点云数量相比较小. 这说明增加的初始姿态的约束造成了全局重建质量的下降. 对比第4、第5列有相同结论. 需要注意的是, 4、5两列虽然视点索引号不同, 但其均是视点优化的结果, 且仅保证了初始位姿是不同的. 比较2、4两列可看出, 考虑权重的情况下, 不论是否约束初始姿态, 重建质量均提高.
3.4 扫描路径与重建结果分析
获得优化的测量视点后, 最终按照式 (14)规划机器人的最短扫描路径. 本文采用插值的方法, 在任意两个视点中插入中间过渡点以保证机器人运动过程不会产生干涉, 最优路径通过自适应游走算法进行计算. 由于本文方法不需要进行现场校准, 被测物体任意放置在测量系统前, 初始姿态具有任意性, 通过表2也说明初始姿态对于重建结果的影响. 所以机器人运动轨迹受到初始姿态的影响, 但是在具有相同的初始姿态下, 机器人轨迹相对固定. 机器人轨迹计算方法也直接用于现场实验, 图7展示了利用本文机器人规划算法设计的仿真界面和现场实验. 图7(a)为仿真界面, 图7(b)为现场实验. 利用本文的测量策略, 在实际测量时, 控制末端相机采图, 并估计目标物体位姿. 然后将该位姿作为初始位姿进行后续视点和路径规划. 规划得到最优路径发送至机器人实现自动测量. 图8(a)展示了两种候选匹配点下的扫描路径.
图 7 机器人扫描规划仿真界面和现场实验Fig. 7 Simulation interface and field experiment of the robot scanning planning改变候选视点数量, 进一步探究不同候选点下扫描路径与重建质量, 探究最优扫描路径在不同候选视点的提升程度. 图8展示了不同候选点集下获得的机器人路径. 不同候选点获得的优化的视点数量不同, 扫描轨迹也不同.
图 8 不同候选匹配点下扫描路径对比Fig. 8 Comparison of scanning paths with different candidate view point由图8可看出, 更多的优化视点获得三维点云数量较多, 重建质量较高, 机器人运动的路径相对较长. 不同的候选点将产生不同的扫描路径, 但本文全部算法流程可根据实际应用中不同初始姿态、不同候选点等情况生成最优的扫描路径. 本文为实现快速摄影测量重建, 重建的点云密度略低于精密摄影测量. 同时, 通过增加候选点数量和重叠点数量(式(11))实现重建质量提高.
利用估计的姿态进行摄影重建, 21个视点重建的时间为90 s. 没有估计姿态直接进行三维重建的时间为111 s, 相对提高了18.92%, 大幅提升了重建速度. 因为没有输入姿态的情况下, 摄影测量系统必须根据特征点估计出相机姿态, 然后进行优化和重建, 故重建速度低. 另外, 结合快速位姿估计、视点优化以及最短路径规划等综合提升了测量系统速度.
3.5 典型特征的摄影测量验证
上文以“金字塔”状凸台零件验证了所提算法. 为验证模型的泛化能力和适用性, 增加球体、柱体、凹台这三类零件的测量位姿估计与视点规划实验. 3D打印件如图9的第1行所示.
图 9 具有典型特征的零件摄影测量位姿估计与视点规划Fig. 9 The pose estimation and viewpoint planning of part with typical features for photogrammetry图9中第1行 ~ 第5行分别展示了3D 打印件的原始图像、位姿估计可视化结果、网络训练结果、最优扫描路径和摄影重建结果. 图9(a) ~ 9(c)表明所提的机器人摄影测量策略适用于常规的零件, 包括球体、柱体和凹台, 具有一定的通用性, 更多的目标物体只需要按照本文方法进行位姿估计、可见性计算、测量规划等工作. 由图中第2行可见可视化的位姿与被测物体十分吻合, 质心坐标重合, 结合第3行测试结果可得出目标包围盒的估计准确率为100%, 同时平移估计精度均优于5 mm, 旋转估计精度优于2°. 第4行给出了最优扫描路径, 其中模型的可见性均满足无重叠投影的要求. 第5行为利用规划的视点进行机器人路径规划得到的摄影重建结果.
同样, 利用估计的位姿进行摄影测量重建, 重建效率对比如表3所示. 本文摄影重建均采用增量式稀疏重建方法, 可缩短整体重建时间, 由于表2验证了熵权对于重建质量的提升, 此处均采用熵权法进行视点优化求解.
表 3 利用深度学习位姿估计的摄影测量效率对比Table 3 Comparison on effectiveness of photogrammetry with estimated pose using deep learning点个数 重建时间 (s) 有初始位姿约束 无初始位姿约束 数量变化 (%) 无位姿估计 有位姿估计 效率提升 (%) 球体 15635 16489 5.18 157 133 15.29 柱体 10138 11503 11.87 182 155 14.84 凹台 11472 12640 9.24 102 83 18.63 由表3可以看出, 考虑初始位姿约束的情况下, 三种模型重建的点云数量均略有减少, 这说明本文方法更加适合稀疏点云的快速重建, 而不能获得稠密的摄影重建. 而融合了位姿估计的摄影测量系统效率均有所提升, 球体、柱体、凹台零件重建效率分别提升15.29%、14.84%和18.63%, 因为利用深度学习的方法进行位姿估计节省了传统摄影测量中位姿依赖本质矩阵分解的计算时间. 综上分析, 融合了初始位姿估计的机器人摄影测量可实现快速稀疏点云重建, 而利用深度学习替代传统依赖特征计算的方法对于摄影测量具有重大意义.
4. 结束语
提出一种基于深度学习初始位姿估计的机器人摄影测量系统视点规划方法, 实现了利用单相机快速位姿估计, 平移精度优于5 mm, 角度精度优于2°. 提出基于相机成像映射的可见性矩阵计算. 利用所提的视点规划策略提升了摄影测量重建的效率. 对几种具有典型特征的零件进行摄影测量, 本文方法在位姿估计、重建效率方面均表现优异, 重建质量因初始位姿约束有所下降, 故适用于快速稀疏摄影重建. 本文方法促进了机器人摄影测量系统现场测量的智能性、灵活性和自动化程度, 对工业摄影测量具有重要的推进意义.
本文网络为开源网络, 为扩展网络应用价值, 公开本文中圆柱体模型的数据集供同行使用. 本文所使用的模型均来自英国诺丁汉大学制造测量团队.
-
图 7 机器人扫描规划仿真界面和现场实验
Fig. 7 Simulation interface and field experiment of the robot scanning planning
图 8 不同候选匹配点下扫描路径对比
Fig. 8 Comparison of scanning paths with different candidate view point
图 9 具有典型特征的零件摄影测量位姿估计与视点规划
Fig. 9 The pose estimation and viewpoint planning of part with typical features for photogrammetry
表 1 信息熵权有效性验证表
Table 1 Effectiveness test for entropy weight
对比实验 目标函数 优化的视
点数(个)点云数(个) A $x^*=\min \displaystyle \sum\limits_{i=1}^N { {w_i}{x_i} }$ 21 21509 B $x^*=\min \displaystyle \sum\limits_{i=1}^N { {x_i} }$ 20 18360 C $x^*=\min \displaystyle \sum\limits_{i=1}^N { {x_i} }$ 21 = 20 (B) +
1 (288)15344 
下载: 导出CSV
表 2 综合权重和初始姿态下重建质量对比
Table 2 Comparison of reconstruction quality with weight and first-sight pose
有权重 无权重 有初始位姿约束 无初始位姿约束 有初始位姿约束 无初始位姿约束 视点索引 1, 13, 100, 113, 143,
149, 173, 189, 190,
196, 207, 269, 272, 28013, 100, 113, 143,
149, 173, 189, 190, 196,
207, 269, 272, 2801, 17, 28, 38, 45, 61,
66, 74, 89, 91, 92, 107,
113, 127, 185, 189,
207, 249, 269, 28014, 35, 43, 45, 56, 59,
73, 75, 89, 111, 127, 149,
162, 185, 189, 207, 249,
256, 274, 281三维点云 点数量 10584 11451 7703 9571 注: 有初始位姿约束下，索引为1的视点需被约束保留.
下载: 导出CSV
表 3 利用深度学习位姿估计的摄影测量效率对比
Table 3 Comparison on effectiveness of photogrammetry with estimated pose using deep learning
点个数 重建时间 (s) 有初始位姿约束 无初始位姿约束 数量变化 (%) 无位姿估计 有位姿估计 效率提升 (%) 球体 15635 16489 5.18 157 133 15.29 柱体 10138 11503 11.87 182 155 14.84 凹台 11472 12640 9.24 102 83 18.63
下载: 导出CSV
-
[1] 郑太雄, 黄帅, 李永福, 冯明驰. 基于视觉的三维重建关键技术研究综述. 自动化学报, 2020, 46(4): 631-652 doi: 10.16383/j.aas.2017.c170502Zheng Tai-Xiong, Huang Shuai, Li Yong-Fu, Feng Ming-Chi. Key techniques for vision based 3D reconstruction: A review. Acta Automatica Sinica, 2020, 46(4): 631-652 doi: 10.16383/j.aas.2017.c170502 [2] 李杰, 李响, 许元铭, 杨绍杰, 孙可意. 工业人工智能及应用研究现状及展望. 自动化学报, 2020, 46(10): 2031-2044 doi: 10.16383/j.aas.200501Lee J, Li Xiang, Xu Yuan-Ming, Yang Shao-Jie, Sun Ke-Yi. Recent advances and prospects in industrial AI and applications. Acta Automatica Sinica, 2020, 46(10): 2031-2044 doi: 10.16383/j.aas.200501 [3] 柴天佑. 工业人工智能发展方向. 自动化学报, 2020, 46(10): 2005-2012 doi: 10.16383/j.aas.c200796Chai Tian-You. Development directions of industrial artificial intelligence. Acta Automatica Sinica, 2020, 46(10): 2005-2012 doi: 10.16383/j.aas.c200796 [4] Sims-Waterhouse D, Bointon P, Piano S, Leach R K. Experimental comparison of photogrammetry for additive manufactured parts with and without laser speckle projection. In: Proceedings of SPIE 10329, Optical Measurement Systems for Industrial Inspection X. Munich, Germany: SPIE, 2017. Article No. 103290W [5] 许杰, 蒋山平, 杨林华, 肖大舟, 张景川. 卫星结构件常压热变形的数字摄影测量. 光学 精密工程, 2012, 20(12): 2667-2673 doi: 10.3788/OPE.20122012.2667Xu Jie, Jiang Shan-Ping, Yang Lin-Hua, Xiao Da-Zhou, Zhang Jing-Shang. Digital photogrammetry for thermal deformation of satellite structures in normal environment. Optics and Precision Engineering, 2012, 20(12): 2667-2673 doi: 10.3788/OPE.20122012.2667 [6] Filion A, Joubair A, Tahan A S, Bonev I A. Robot calibration using a portable photogrammetry system. Robotics and Computer-Integrated Manufacturing, 2018, 49: 77-87 doi: 10.1016/j.rcim.2017.05.004 [7] Kinnell P, Rymer T, Hodgson J, Justham L, Jackson M. Autonomous metrology for robot mounted 3D vision systems. CIRP Annals, 2017, 66(1): 483-486 doi: 10.1016/j.cirp.2017.04.069 [8] Kwon H, Na M, Song J B. Rescan strategy for time efficient view and path planning in automated inspection system. International Journal of Precision Engineering and Manufacturing, 2019, 20(10): 1747-1756 doi: 10.1007/s12541-019-00186-x [9] Raffaeli R, Mengoni M, Germani M, Mandorli F. Off-line view planning for the inspection of mechanical parts. International Journal on Interactive Design and Manufacturing (IJIDeM), 2013, 7(1): 1-12 doi: 10.1007/s12008-012-0160-1 [10] Li L N, Xu D, Niu L K, Lan Y, Xiong X Y. A path planning method for a surface inspection system based on two-dimensional laser profile scanner. International Journal of Advanced Robotic Systems, 2019, 16(4): Article No. 1729881419862463 [11] Alsadik B, Gerke M, Vosselman G. Visibility analysis of point cloud in close range photogrammetry. In: Proceedings of the ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences. Riva del Garda, Italy: ISPRS, 2014. 9−16 [12] Tarbox G H, Gottschlich S N. Planning for complete sensor coverage in inspection. Computer Vision and Image Understanding, 1995, 61(1): 84-111 doi: 10.1006/cviu.1995.1007 [13] Jing W, Polden J, Lin W, Shimada K. Sampling-based view planning for 3D visual coverage task with unmanned aerial vehicle. In: Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). Daejeon, South Korea: IEEE, 2016. 1808−1815 [14] Jing W, Polden J, Tao P Y, Goh C F, Lin W, Shimada K. Model-based coverage motion planning for industrial 3D shape inspection applications. In: Proceedings of the 13th IEEE Conference on Automation Science and Engineering (CASE). Xi＇an, China: IEEE, 2017. 1293−1300 [15] 姜涛, 程筱胜, 崔海华, 田威. 面向机器人位姿测量的大视场变焦测量方法. 光学学报, 2018, 38(8): Article No. 0815012Jiang Tao, Cheng Xiao-Sheng, Cui Hai-Hua, Tian Wei. Large field of view vision method for robot pose measurement based on zoom lens. Acta Optica Sinica, 2018, 38(8): Article No. 0815012 [16] Jiang T, Cheng X, Cui H, Li X. Combined shape measurement based on locating and tracking of an optical scanner. Journal of Instrumentation, 2019, 14(1): Article No. P01006 [17] Tekin B, Sinha S N, Fua P. Real-time seamless single shot 6D object pose prediction. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. Salt Lake City, USA: IEEE, 2018. 292−301 [18] Su H, Qi C R, Li Y Y, Guibas L J. Render for CNN: Viewpoint estimation in images using CNNs trained with rendered 3D model views. In: Proceedings of the IEEE International Conference on Computer Vision (ICCV). Santiago, Chile: IEEE, 2015. 2686−2694 [19] Tulsiani S, Malik J. Viewpoints and keypoints. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR). Boston, USA: IEEE, 2015. 1510−1519 [20] Kendall A, Grimes M, Cipolla R. PoseNet: A convolutional network for real-time 6-DOF camera relocalization. In: Proceedings of the IEEE International Conference on Computer Vision (ICCV). Santiago, Chile: IEEE, 2015. 2938−2946 [21] Xiang Y, Schmidt T, Narayanan V, Fox D. PoseCNN: A convolutional neural network for 6D object pose estimation in cluttered scenes. arXiv preprint arXiv: 1711.00199, 2017. [22] Kehl W, Manhardt F, Tombari F, Ilic S, Navab N. SSD-6D: Making RGB-based 3D detection and 6D pose estimation great again. In: Proceedings of the IEEE International Conference on Computer Vision (ICCV). Venice, Italy: IEEE, 2017. 1530−1538 [23] Rad M, Lepetit V. BB8: A scalable, accurate, robust to partial occlusion method for predicting the 3D poses of challenging objects without using depth. In: Proceedings of the IEEE International Conference on Computer Vision (ICCV). Venice, Italy: IEEE, 2017. 3848−3856 [24] Carpin S, Pillonetto G. Motion planning using adaptive random walks. IEEE Transactions on Robotics, 2005, 21(1): 129-136 doi: 10.1109/TRO.2004.833790 [25] Cortsen J, Petersen H G. Advanced off-line simulation framework with deformation compensation for high speed machining with robot manipulators. In: Proceedings of the IEEE/ASME International Conference on Advanced Intelligent Mechatronics (AIM). Kaohsiung, China: IEEE, 2012. 934−939 期刊类型引用(0)
其他类型引用(3)
-
下载:
计量
- 文章访问数: 1577
- HTML全文浏览量: 1036
- PDF下载量: 341
- 被引次数: 3
