Investigation on mechanism of material ejection by nanosecond laser ablation of mono-crystalline silicon under different environments
-
摘要: 为了探究纳秒激光烧蚀单晶硅过程中熔融物质的喷溅过程, 对不同加工环境中纳秒激光烧蚀单晶硅的过程进行了模拟仿真及实验验证。采用Level-Set界面追踪法, 通过仿真软件建立有限元模型, 对空气、静水和真空加工环境中纳秒激光烧蚀单晶硅时物质抛出过程进行模拟仿真, 研究了不同加工环境中温度场、速度场对表面喷溅的影响; 采用波长266 nm、脉宽30 ns和频率50 Hz的单脉冲激光烧蚀单晶硅的工艺实验对仿真结果进行验证, 通过原子力显微镜和数字显微镜对烧蚀结果进行了表征。结果表明, 空气环境中, t=30 ns时熔融物质的喷溅速率达到14.1 m/s, 在微孔内部蒸汽压力的作用下, 烧蚀区域熔融物质向外喷出; 静水环境中, t=30 ns时熔融物质的喷溅速率为1.68 m/s, 远低于空气中的喷溅速率, 熔融物质快速冷却; 真空环境中, 材料在短时间内汽化, t=30 ns时熔融物质的喷溅速率最大可达18.4 m/s, 较高的喷溅速率有利于物质的抛出。加工环境对纳秒激光烧蚀单晶硅的物质抛出影响较大, 这为提高纳秒激光加工单晶硅的加工质量提供了参考。Abstract: In order to investigate the process of molten material ejection during nanosecond laser ablation of mono-crystalline silicon, simulation and experimental investigation on nanosecond laser ablation of silicon under different environments were carried out. Simulation on the material ejection during nanosecond laser ablation of silicon under air, water, and vacuum environments with simulation software was carried out based on the Level-Set method. The effects of the temperature field and velocity field on surface sputtering under different environments were investigated. Experiments on laser ablation of mono-crystalline silicon using a wavelength of 266 nm, a pulse width of 30 ns, and a frequency of 50 Hz were done to verify the simulation results. The ablation results were observed by atomic force microscopy and digital microscopy. The results show that the spatter velocity under the air environment is 14.1 m/s at t=30 ns, and the molten material ejectes the ablated holes outwards due to the action of vapor pressure in the ablated area. Under the water environment, the spatter velocity is 1.68 m/s at t=30 ns, much lower than that in air conditions. Under the vacuum environment, the material is vaporized for a short time, and the spatter velocity of the molten material is 18.4 m/s at t=30 ns. The molten materials eject productively due to the high spatter velocity. The adhesion of molten material on the surface is improved. The processing environment has a strong influence on the material ejection during the nanosecond laser ablation of silicon. The investigation provides a reference for improving the quality during nanosecond laser micro-machining.
-
0. 引言
利用图像全站仪技术进行自动工业摄影测量在技术上已成为可能。采用图像全站仪进行摄影测量具备高精度、无接触、面测量、自动化及无人为照准误差等特性，在天文测量、变形测量及精密工程测量等领域有广泛应用前景[1-7]。高精度的相机标定是图像全站仪进行精密摄影测量一个基础而又不可或缺的问题，也是目前图像全站仪及其应用亟待解决的关键技术之一，但目前相关的研究鲜见报道，相近的研究也主要集中在天文测量和变形测量领域，且只针对望远镜相机定焦测量展开研究，未对广角相机展开研究[8-14]。
EHRHART等人在使用徕卡MS50图像全站仪进行工程结构动态监测时，对该望远镜相机检校进行了研究，认为对于变形监测而言，可以不考虑相机的像主点坐标和畸变，只需要相机的角分辨率即可[11]。国内信息工程大学天文课题组在使用TS50i图像全站仪进行天文测量时，对相机的像平面倾斜、像主点坐标、主距及畸变系数等进行了检校[12]，提出的检校方法简单且方便操作，但相机成像模型并不严密，未考虑相机在全站仪设备中的偏心元素(包括3个位置偏移量x0, y0, z0和3个角度偏移量α, β, γ)，其本质是将全站仪坐标系统等价于相机坐标系统，这会引起相机标定误差。对于望远镜相机而言，相机与全站仪望远镜同轴安置，位置偏移量级为亚毫米级，角度偏移量级为角秒级，全站仪仪器中心近似为相机光心对相机标定结果影响较小，但对于广角相机而言，相机与全站仪望远镜显著偏心安置，仪器中心与相机光心的位置偏差为厘米级，将仪器中心当作相机光心会显著引起相机标定误差，使得相机标定参数在不同摄影距离下标定解算的结果不具备一致性，无法实现图像测量和全站仪角度测量的精确统一，从而限制了图像全站仪在精密工程测量领域中的应用。综上所述，对图像全站仪相机科学的标定方法展开研究具有应用意义。
偏心元素的求解是图像全站仪相机标定独有的问题，由于图像全站仪相机与全站仪刚性连接，偏心元素的各项数值均为常数，在一般标定方法中存在参数耦合问题，需用特殊的方法进行参数解耦。ZHANG提出使用检校条，通过多测站近距离旋转摄影的方式，构建成像共线方程，采用摄影测量前方交会的方式完成偏心元素检校[15]。ZHU提出相机参数逐级标定法，利用全站仪做小角旋转单独标定相机主距，将之作为已知量逐步求解其它参数，避免参数之间的耦合[16]。两者均从理论上指出图像全站仪可通过旋转摄影的方式构建标定场，但缺乏实验来验证旋转摄影标定方法的有效性。
本文作者基于旋转摄影理论，采用图像全站仪单目标旋转摄影标定的方法，推导并建立了图像全站仪相机标定模型，同时，开展实验对模型在不同摄影距离上进行了相机标定，均获得了高精度的标定结果。
1. 单目标旋转摄影标定方法
本文中提出的单目标旋转摄影标定方法主要是指首先通过计算机程序控制全站仪相机对单个靶标进行扫描式摄影，然后利用相机摄影时的姿态角及对应靶标影像的中心像素坐标来构建相机成像方程，最后利用构建的成像方程实现对相机标定参数的求解。该方法不仅充分发挥了图像全站仪平台特性，同时对目标对象的布设以及标定过程的自动化程度方面，相对于传统的方法，具备无可比拟的优势。
单目标旋转摄影标定方法中，人们借助图像全站仪的GEOCOM二次开发接口[17]，通过编写计算机程序实现靶标图像自动化数据采集，整个测量和拍摄作业自动化程度高。具体流程是：通过计算机控制指令驱动望远镜在水平和垂直方向上指向规划的角度并记录当前望远镜姿态，相机与望远镜固联，随全站仪旋转驱动；摄影标定目标不动，仪器每次旋转规划的角度后，对靶标进行拍摄，使目标近似均匀的成像在电荷耦合器件(charge-coupled device, CCD)不同区域。该过程可等价成相机不动而靶标在以相机为中心的球面上移动[18]，目标的每一个姿态角均可等价于一个方向控制点，且控制点精度仅取决于测角精度。在几十米范围内，图像全站仪测角精度可达到0.5″，远高于测距精度，这使得旋转摄影标定方法获得的精度比控制网标定方法更高。采用单靶标进行拍照，可将靶标做得相对较大，有利于提升靶标中心提取精度。用旋转摄影的方法进行点位扩展，相当于在空间中布设高精度的异面控制点，能近乎零成本和更高效地建立高精度且视域内均匀分布的空间标定场。
2. 单目标旋转摄影标定的理论模型
2.1 坐标系
相机标定过程中，涉及到不同坐标系间的坐标转换，主要有全站仪坐标系、望远镜坐标系、相机坐标系以及像素坐标系。为保持与全站仪观测量的一致性，均采用了与全站仪坐标系一致的左手坐标系。
(a) 全站仪坐标系(total station，TS)OTS-XTSYTSZTS，原点为全站仪中心，ZTS轴取全站仪竖轴，XTS轴平行水平度盘零方向，YTS轴由左手系确定。
(b) 望远镜坐标系(telescope，T)OT-XTYTZT，原点取全站仪中心，ZT轴取全站仪望远镜视准轴，YT轴取全站仪横轴，XT轴由左手系确定。当望远镜坐标系处于置零状态时与全站仪坐标系一致。
(c) 相机坐标系(camera，C)OC-XCYCZC，原点为相机光心，ZC轴取相机主光轴，XC轴与像素坐标系XP轴平行且同方向，YC轴由左手系确定。
(d) 像素坐标系(pixel，P)OP-XPYP，建立在相机坐标系ZC=-f的2维平面坐标系，f为相机等效主距。原点为图像左上角，取行方向为XP轴，列方向为YP轴。
图像全站仪相机与全站仪刚性连接，偏心元素为常量，相机坐标系与望远镜坐标系的空间关系如图 1所示。各坐标系的空间关系示意图如图 2所示，其中相机坐标系的ZC轴延长线与像素坐标系平面的交点为像主点(u0, v0)。
![图 1 T坐标系与C坐标系空间关系图]()
图 1 T坐标系与C坐标系空间关系图Figure 1. Spatial relationship between T-coordinate system and C-coordinate system2.2 相机成像模型
首先，建立T坐标系至S坐标系的位姿转换。设全站仪照准靶标P时，全站仪水平角、天顶角和距离读数分别为θ0, φ0, s，P点在全站仪坐标系下的坐标为：
\boldsymbol{P}_{\mathrm{TS}}=s\left[\begin{array}{c} \sin \varphi_{0} \cos \theta_{0} \\ \sin \varphi_{0} \sin \theta_{0} \\ \cos \varphi_{0} \end{array}\right] (1) 驱动全站仪旋转到某一规划指向，此时全站仪水平角、天顶角读数为分别为θ, φ。根据TS坐标系与T坐标系的空间旋转关系可知，TS坐标系绕ZTS轴旋转θ，再绕YTS轴旋转φ，可实现TS坐标系到T坐标系的转换，即：
\begin{equation*} \boldsymbol{P}_{\mathrm{T}}=\boldsymbol{R}_{Y_{\mathrm{TS}}}(\varphi) \cdot \boldsymbol{R}_{Z_{\mathrm{TS}}}(\theta) \cdot \boldsymbol{X}_{\mathrm{TS}} \end{equation*} (2) 式中: RZTS(θ)、RYTS(φ)为对应旋转矩阵，可表示为：
\begin{align*} & \boldsymbol{R}_{Z_{\mathrm{TS}}}(\theta)=\left[\begin{array}{ccc} \cos \theta & \sin \theta & 0 \\ -\sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{array}\right] \end{align*} (3) \begin{align*} & \boldsymbol{R}_{Y_{\mathrm{TS}}}(\varphi)=\left[\begin{array}{ccc} \cos \varphi & 0 & -\sin \varphi \\ 0 & 1 & 0 \\ \sin \varphi & 0 & \cos \varphi \end{array}\right] \end{align*} (4) 式(2)可表示为：
\begin{equation*} \boldsymbol{P}_{\mathrm{T}}=s \cdot \boldsymbol{n}_{\mathrm{T}} \end{equation*} (5) 式中: nT为靶标中心在望远镜坐标系下的单位方向向量，其计算表达式为：
\begin{gather*} \boldsymbol{n}_{\mathrm{T}}=\left[\begin{array}{l} s_{x} \\ s_{y} \\ s_{z} \end{array}\right]= \\ {\left[\begin{array}{c} \cos \varphi \sin \varphi_{0} \cos \left(\theta_{0}-\theta\right)-\sin \varphi \cos \varphi_{0} \\ \sin \varphi_{0} \sin \left(\theta_{0}-\theta\right) \\ \sin \varphi \sin \varphi_{0} \cos \left(\theta_{0}-\theta\right)+\cos \varphi \cos \varphi_{0} \end{array}\right]} \end{gather*} (6) 其次，建立T坐标系至C坐标系的位姿转换，根据T坐标系与C坐标系的空间关系有：
\boldsymbol{P}_{\mathrm{C}}=\boldsymbol{R}_{X_{\mathrm{T}}}(\alpha) \cdot \boldsymbol{R}_{Y_{\mathrm{T}}}(\beta) \cdot \boldsymbol{R}_{Z_{\mathrm{T}}}(\gamma) \cdot \boldsymbol{P}_{\mathrm{T}}+\left[\begin{array}{l} x_{0} \\ y_{0} \\ z_{0} \end{array}\right] $$ (7) 式中: α, β, γ为T坐标系与C坐标系的角度偏移量; x0, y0, z0为位置偏移量。
考虑到无论是广角相机或是望远镜相机，C坐标系至T坐标系的空间转换均为小角度旋转，可将旋转矩阵线性化并取一次项近似代替，有：
\boldsymbol{P}_{\mathrm{C}}=\left[\begin{array}{ccc} 1 & \gamma & -\beta \\ -\gamma & 1 & \alpha \\ \beta & -\alpha & 1 \end{array}\right] \boldsymbol{P}_{\mathrm{T}}+\left[\begin{array}{l} x_{0} \\ y_{0} \\ z_{0} \end{array}\right] (8) 式(8)也可写为：
\begin{equation*} \boldsymbol{P}_{\mathrm{C}}=s \cdot \boldsymbol{n}_{\mathrm{C}} \end{equation*} (9) 式中: nC为靶标中心在相机坐标系下的单位方向向量, 其计算表达式为：
\boldsymbol{n}_{\mathrm{C}}=\left[\begin{array}{c} c_{x} \\ c_{y} \\ c_{z} \end{array}\right]=\left[\begin{array}{c} s_{x}+\gamma s_{y}-\beta s_{z}+x_{0} / s \\ -\gamma s_{x}+s_{y}+\alpha s_{z}+y_{0} / s \\ \beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s \end{array}\right] (10) 最后，建立C坐标系至P坐标系的坐标转换，由相机成像共线方程有：
\begin{equation*} \frac{x-u_{0}}{c_{x}}=\frac{y-v_{0}}{c_{y}}=\frac{f}{c_{z}} \end{equation*} (11) 由式(11)可得图像全站仪成像模型为：
\left\{\begin{array}{l} x=f_{x} \frac{c_{x}}{c_{z}}+u_{0} \\ y=f_{y} \frac{c_{y}}{c_{z}}+v_{0} \end{array}\right. (12) 式中: x, y为靶标中心像素坐标; fx, fy为相机等效主距; u0, v0为像主点坐标。式(12)也可表示为：
\left\{\begin{array}{l} x=f_{x} \frac{s_{x}+\gamma s_{y}-\beta s_{z}+x_{0} / s}{\beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s}+u_{0} \\ y=f_{y} \frac{-\gamma s_{x}+s_{y}+\alpha s_{z}+y_{0} / s}{\beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s}+v_{0} \end{array}\right. (13) 2.3 相机标定方程
对式(13)进行线性化，可得线性化观测误差方程为：
\begin{equation*} \boldsymbol{V}=\boldsymbol{A} \boldsymbol{X}-\boldsymbol{L} \end{equation*} (14) 式中: V为观测残差矩阵; A为参数系数矩阵; X为待解的相机参数; L为观测矩阵。其中：
\begin{gather*} \boldsymbol{X}=\left[u_{0}, v_{0}, f_{x}, f_{y}, \alpha, \beta, \gamma, x_{0}, y_{0}, z_{0}\right]^{\mathrm{T}} \end{gather*} (15) \boldsymbol{L}=\left[\begin{array}{l} x \\ y \end{array}\right] (16) 相机参数系数矩阵的计算表达式为：
\boldsymbol{A}=\left[\begin{array}{l} \frac{\partial x}{\partial u_{0}} \frac{\partial x}{\partial v_{0}} \frac{\partial x}{\partial f_{x}} \frac{\partial x}{\partial f_{y}} \frac{\partial x}{\partial \alpha} \frac{\partial x}{\partial \beta} \frac{\partial x}{\partial \gamma} \frac{\partial x}{\partial x_{0}} \frac{\partial x}{\partial y_{0}} \frac{\partial x}{\partial z_{0}} \\ \frac{\partial y}{\partial u_{0}} \frac{\partial y}{\partial v_{0}} \frac{\partial y}{\partial f_{x}} \frac{\partial y}{\partial f_{y}} \frac{\partial y}{\partial \alpha} \frac{\partial y}{\partial \beta} \frac{\partial y}{\partial \gamma} \frac{\partial y}{\partial x_{0}} \frac{\partial y}{\partial y_{0}} \frac{\partial y}{\partial z_{0}} \end{array}\right] (17) 对式(17)求导取1阶项，简化整理可得：
\mathit{\boldsymbol{A = }}\left[ \begin{array}{l} 1\quad 0\quad\frac{{{c_x}}}{{{c_z}}}\quad 0\quad \quad \left( {x - {u_0}} \right)\frac{{{s_y}}}{{{c_z}}}\quad \quad \left( {{u_0} - x} \right)\left( {\frac{{{s_x}}}{{{c_z}}} + \frac{{{s_z}}}{{{c_x}}}} \right)\quad \left( {x - {u_0}} \right)\frac{{{s_y}}}{{{c_x}}}\quad \frac{{x - {u_0}}}{{s{c_x}}}\quad 0\quad \frac{{{u_0} - x}}{{s{c_z}}}\\ 0\quad 1\quad 0\quad \frac{{{c_y}}}{{{c_z}}}\quad \left( {{v_0} - y} \right)\left( {\frac{{{s_z}}}{{{c_y}}} + \frac{{{s_y}}}{{{c_z}}}} \right)\quad \quad \left( {y - {v_0}} \right)\frac{{{s_x}}}{{{c_z}}}\quad\quad \left( {{v_0} - y} \right)\frac{{{s_x}}}{{{c_y}}}\quad \quad 0\quad \frac{{y - {v_0}}}{{s{c_y}}}\quad \frac{{{v_0} - y}}{{s{c_z}}} \end{array} \right] (18) 采用最小二乘方法求解相机参数，有：
\begin{equation*} \boldsymbol{X}=\left(\boldsymbol{A}^{\mathrm{T}} \boldsymbol{P} \boldsymbol{A}\right)^{-1} \boldsymbol{A} \boldsymbol{P} \boldsymbol{L} \end{equation*} (19) 式中: P为观测量的权矩阵。由于标定参数间存在的强相关性，直接组建标定方程求解时标定方程可能奇异，无法正确求解，应当分析变量间的相关性，解决参数强相关性问题。
2.4 相关性分析
通过参数之间的偏导数关系，对参数的相关性展开分析，具体为：
\begin{gather*} {\left[\begin{array}{llll} \frac{\partial \beta}{\partial x_{0}} & \frac{\partial \beta}{\partial u_{0}} & \frac{\partial x_{0}}{\partial u_{0}} & \frac{\partial z_{0}}{\partial f_{x}} \\ \frac{\partial \alpha}{\partial y_{0}} & \frac{\partial \alpha}{\partial v_{0}} & \frac{\partial y_{0}}{\partial v_{0}} & \frac{\partial z_{0}}{\partial f_{y}} \end{array}\right]=} \\ {\left[\begin{array}{llll} \frac{\partial x / \partial x_{0}}{\partial x / \partial \beta} & \frac{\partial x / \partial u_{0}}{\partial x / \partial \beta} & \frac{\partial x / \partial u_{0}}{\partial x / \partial x_{0}} & \frac{\partial x / \partial f_{x}}{\partial x / \partial z_{0}} \\ \frac{\partial y / \partial y_{0}}{\partial y / \partial \alpha} & \frac{\partial y / \partial v_{0}}{\partial y / \partial \alpha} & \frac{\partial y / \partial v_{0}}{\partial y / \partial y_{0}} & \frac{\partial y / \partial f_{y}}{\partial y / \partial z_{0}} \end{array}\right]} \end{gather*} (20) 经计算整理有：
\begin{array}{c} {\left[\begin{array}{llll} \frac{\partial \beta}{\partial x_{0}} & \frac{\partial \beta}{\partial u_{0}} & \frac{\partial x_{0}}{\partial u_{0}} & \frac{\partial z_{0}}{\partial f_{x}} \\ \frac{\partial \alpha}{\partial y_{0}} & \frac{\partial \alpha}{\partial v_{0}} & \frac{\partial y_{0}}{\partial v_{0}} & \frac{\partial z_{0}}{\partial f_{y}} \end{array}\right]=} \\ {\left[\begin{array}{cccc} \frac{-c_{z} / s}{c_{x} s_{x}+c_{z} s_{z}} & \frac{-c_{z}^{2} / f_{x}}{c_{x} s_{x}+c_{z} s_{z}} & \frac{s c_{z}}{f_{x}} & \frac{s c_{z}}{-f_{x}} \\ \frac{c_{z} / s}{c_{y} s_{y}+c_{z} s_{z}} & \frac{c_{z}^{2} / f_{y}}{c_{y} s_{y}+c_{z} s_{z}} & \frac{s c_{z}}{f_{y}} & \frac{s c_{z}}{-f_{y}} \end{array}\right]} \end{array} (21) 式中: 相机主距fx和fy可视作常数，偏心元素的各项数值均为小量常数，可计算出sz, cz的值趋于1，远大于sx, sy, cx, cy。当考虑单距离旋转摄影标定时，即摄影距离s为常量时，式(21)中各项导数计算值均近似为常数，即β, x0, u0强相关，α、y0、v0强相关，fx, fy, z0强相关，故通过单距离观测进行相机标定时，仅可考虑变量u0, v0, fx, fy, γ进行标定解算，无法解算出参数α, β, x0, y0, z0。
当考虑多距离旋转摄影标定时，即摄影距离s为变量时，式(21)中∂β/∂u0, ∂α/∂v0的值近似不变，其余项的计算值线性变化，即x0, y0, z0与摄影距离s线性相关，故通过多距离观测整体解算可解耦x0, y0, z0，但无法解算α, β。经分析认为，需改变相机相对全站仪的位置才能解耦α, β，这可以通过切换全站仪度盘来改变相机相对全站仪的位置，采用双盘观测的方式解耦α, β。
2.5 双盘观测方程
将望远镜绕水平方向旋转180°，再绕竖直方向旋转180°可实现全站仪度盘切换。切换观测度盘会改变相机相对全站仪的位置，假设在盘左观测中偏心元素为x0, y0, z0, α, β, γ，在盘右观测中偏心元素则为-x0, -y0, z0, -α, -β, γ。令盘左观测时靶标像素坐标为(xl, yl)，盘右观测时靶标像素坐标为(xr, yr)。盘右图像成正像显示，观测数据需作镜像处理，将式(12)对应参数进行修正，可得盘右观测时对应的成像方程为：
\left\{\begin{array}{l} x_{\mathrm{r}}=f_{x} \frac{-s_{x}-\gamma s_{y}-\beta s_{z}+x_{0} / s}{-\beta s_{x}+\alpha s_{y}+s_{z}+z_{0} / s}+\left(R-u_{0}\right) \\ y_{\mathrm{r}}=f_{y} \frac{\gamma s_{x}-s_{y}+\alpha s_{z}+y_{0} / s}{-\beta s_{x}+\alpha s_{y}+s_{z}+z_{0} / s}+\left(C-v_{0}\right) \end{array}\right. (22) 式中: R为图像的行数; C为图像的列数。构建双盘观测方程为：
\left\{\begin{array}{l} x_{1}=f_{x} \frac{s_{x}+\gamma s_{y}-\beta s_{z}+x_{0} / s}{\beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s}+u_{0} \\ y_{1}=f_{y} \frac{-\gamma s_{x}+s_{y}+\alpha s_{z}+y_{0} / s}{\beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s}+v_{0} \\ x_{\mathrm{r}}=f_{x} \frac{-s_{x}-\gamma s_{y}-\beta s_{z}+x_{0} / s}{-\beta s_{x}+\alpha s_{y}+s_{z}+z_{0} / s}+\left(R-u_{0}\right) \\ y_{\mathrm{r}}=f_{y} \frac{\gamma s_{x}-s_{y}+\alpha s_{z}+y_{0} / s}{-\beta s_{x}+\alpha s_{y}+s_{z}+z_{0} / s}+\left(C-v_{0}\right) \end{array}\right. (23) 利用式(23)，在多个距离上扫描目标的盘左和盘右成像坐标，通过整体平差解算可求解出全部相机参数。借鉴分步求解思路来降低参数耦合的影响，先假设参数α、β值为0°，通过整体平差解出参数fx, fy, u0, v0, γ, x0, y0, z0，此时α、β耦合在u0、v0中，后通过某种方式单独解算α、β，再对u0、v0进行修正。
由式(23)可构建双盘观测量：
\left\{\begin{align*} \frac{x_{1}+x_{\mathrm{r}}}{2}= & \frac{f_{x}}{2}\left(\frac{s_{x}+\gamma s_{y}-\beta s_{z}+x_{0} / s}{\beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s}+\right. \\ & \left.\frac{-s_{x}-\gamma s_{y}-\beta s_{z}+x_{0} / s}{-\beta s_{x}+\alpha s_{y}+s_{z}+z_{0} / s}\right)+\frac{R}{2} \\ \frac{y_{1}+y_{\mathrm{r}}}{2}= & \frac{f_{y}}{2}\left(\frac{-\gamma s_{x}+s_{y}+\alpha s_{z}+y_{0} / s}{\beta s_{x}-\alpha s_{y}+s_{z}+z_{0} / s}+\right. \\ & \left.\frac{\gamma s_{x}-s_{y}+\alpha s_{z}+y_{0} / s}{-\beta s_{x}+\alpha s_{y}+s_{z}+z_{0} / s}\right)+\frac{C}{2} \end{align*}\right. (24) 忽略小量βsx, αsy，整理化简可得：
\left\{\begin{array}{l} \frac{x_{1}+x_{\mathrm{r}}}{2}=-\beta f_{x}+f_{x} \frac{x_{0}}{s s_{z}}+\frac{R}{2} \\ \frac{y_{1}+y_{\mathrm{r}}}{2}=\alpha f_{y}+f_{y} \frac{y_{0}}{s s_{z}}+\frac{C}{2} \end{array}\right. (25) 可得到α, β的计算式为:
\left\{\begin{array}{l} \beta=\frac{x_{0}}{s s_{z}+z_{0}}+\frac{R}{2 f_{x}}-\frac{x_{1}+x_{\mathrm{r}}}{2 f_{x}} \\ \alpha=-\frac{y_{0}}{s s_{z}+z_{0}}-\frac{C}{2 f_{y}}+\frac{y_{1}+y_{\mathrm{r}}}{2 f_{y}} \end{array}\right. (26) 在解得α, β的基础上，对像主点坐标做相应修正，可获得实际像主点坐标值u0, v0：
\left\{\begin{array}{l} \bar{u}_{0}=u_{0}-f_{x} \beta \\ \bar{v}_{0}=v_{0}+f_{y} \alpha \end{array}\right. (27) 2.6 相机畸变模型
由于镜头设计的复杂性和工艺水平等因素的影响，镜头折射率不一致，光学器件存在偏移，导致实际成像点(x′, y′)与理论像点(x, y)之间存在像素偏差(δx, δy)，即：
\left[\begin{array}{l} x \\ y \end{array}\right]=\left[\begin{array}{l} x^{\prime} \\ y^{\prime} \end{array}\right]+\left[\begin{array}{l} \delta_{x} \\ \delta_{y} \end{array}\right] (28) Brown-Conrady畸变理论认为存在3种畸变类型，分别为径向畸变、偏心畸变和薄棱镜畸变[19]。同时考虑这3种畸变，相机畸变模型可表示为：
\left\{\begin{align*} \delta_{x}= & \left(x-u_{0}\right)\left(k_{1} r^{2}+k_{2} r^{4}+k_{3} r^{6}\right)+ \\ & 2 \rho_{1}\left(x-u_{0}\right)\left(y-v_{0}\right)+\rho_{2}\left(r^{2}+2 x^{2}\right)+s_{1} r^{2} \\ \delta_{y}= & \left(y-v_{0}\right)\left(k_{1} r^{2}+k_{2} r^{4}+k_{3} r^{6}\right)+ \\ & 2 \rho_{2}\left(x-u_{0}\right)\left(y-v_{0}\right)+\rho_{1}\left(r^{2}+2 y^{2}\right)+s_{2} r^{2} \end{align*}\right. (29) 式中: r=\sqrt{\left(x-u_0\right)^2+\left(y-v_0\right)^2}，表示像点到透视中心的像素距离；k1, k2, k3为径向畸变系数；ρ1和ρ2为偏心畸变系数；s1和s2为薄棱镜畸变系数。
3. 标定实验及分析
3.1 实验描述
基于上述理论，本文作者对TS60图像全站仪进行相机标定实验。TS60同时具有望远镜相机和广角相机，望远镜相机为变焦相机，视场角为1.5°，广角相机为定焦相机，视场角为19.4°，成像规格为2560 pixel×1920 pixel(R×C)。标定靶标采用同心圆环靶标，设计有十字丝便于全站仪初始精确照准。设计在8 m、12 m、16 m、20 m、24 m、30 m的距离进行望远镜相机和广角相机的图像采集。设定合适的摄影间隔使得标定靶标在视域内均匀分布成像，本文中设计旋转摄影照片数量为26×24，对应广角相机旋转摄影方位角和高度角摄影间隔分别为35.77′和29.25′，对应望远镜相机旋转摄影方位角和高度角摄影间隔分别为2.77′和2.25′，在单个距离采集624张靶标图像。基于C#实现计算机与全站仪的蓝牙通讯，参照手册开发计算机程序控制全站仪相机驱动到规划的角度姿态对靶标进行拍照，实现靶标图像自动化数据采集。实验室现场环境如图 3所示。望远镜相机获得的部分靶标图像集如图 4所示。
相机标定流程如下：(a)通过中心提取算法[20-21]获取靶标中心像素坐标(x, y)，并记录对应的摄影方向(θ, φ)；(b)仅考虑参数u0, v0, fx, fy, γ进行解算，作为后续求解的初值；(c)采用多距离双盘观测整体解算模型解算变量u0, v0, fx, fy, γ, x0, y0, z0；(d)构建双盘观测量解算α, β，在此基础上对像主点坐标进行改正；(e)广角相机考虑相机畸变参数k1、k2，望远镜相机考虑畸变参数k1，求解相机畸变参数；(f)通过畸变参数对靶标像素坐标做相应改正，重复步骤(b)~(f)，直至求解的参数变化量小于容许值，获得最终相机标定参数。
3.2 实验结果及分析
表 1中为采用单个距离的望远镜相机盘左和盘右观测数据进行相机标定的结果。表中, 全站仪观测度盘面值(face)为0表示盘左观测，度盘面值为1表示盘右观测，mr表示平均重投影点位误差，为中心提取算法获取的靶标中心像素坐标与通过模型计算得到的靶标中心像素坐标的平均差值。表 2中为联合所有观测距离的望远镜相机盘左和盘右观测数据分别进行标定解算的结果。
表 1 望远镜相机单距离观测标定解算Table 1. Calibration solution of telescope camera for single distance observationface s/m u0/pixel v0/pixel fx /pixel fy /pixel γ /(″) mr/pixel 0 8.11 962.47 1280.89 105588.46 105588.45 487.45 0.12 1 8.11 958.52 1279.37 105565.96 105592.29 491.36 0.10 0 12.27 962.05 1281.44 105484.74 105526.00 490.01 0.10 1 12.27 961.19 1280.47 105492.32 105529.98 478.04 0.12 0 15.83 962.02 1281.02 105512.10 105550.65 492.63 0.08 1 15.83 958.84 1278.62 105512.43 105542.27 477.48 0.07 0 19.85 962.37 1283.05 105522.13 105499.94 461.41 0.37 1 19.85 959.45 1279.93 105501.19 105508.14 489.65 0.14 0 23.93 960.35 1283.33 105465.85 105487.35 481.27 0.12 1 23.93 961.68 1280.78 105471.87 105584.78 522.80 0.18 0 29.70 961.57 1283.24 105462.67 105444.87 475.50 0.26 1 29.70 958.84 1277.76 105447.33 105467.91 450.16 0.19 表 2 望远镜相机多距离观测标定解算Table 2. Calibration solution of telescope camera for multiple distance observationface u0/pixel v0/pixel fx/pixel fy/pixel α/(″) β/(″) γ/(″) x0/mm y0/mm z0/mm k1/10-10 mr/pixel 0 960.63 1282.91 105392.33 105382.56 -6.30 2.30 496.28 0.16 -0.15 -15.77 -3.13 0.88 1 959.46 1279.69 105412.70 105452.90 -6.30 2.30 493.64 0.01 -0.01 -12.13 -4.76 1.15 由表 1与表 2可知：望远镜相机不同距离同盘观测的像主点坐标较为稳定，同一距离不同度盘观测的像主点坐标平均差异为2 pixel；单距离观测相机主距有随摄影距离增加而减小的趋势，与多距离观测主距最大差异为205 pixel，对应引起图像边缘的角度最大变化约7.0″；单盘解算平均重投影误差为0.15 pixel，对应图像分辨率0.3″，多距离解算平均重投影误差为1.02 pixel，对应的图像角分辨率约2.0″。
表 3中为采用单个距离的广角相机盘左和盘右观测数据进行相机标定的结果。表 4中为联合所有观测距离的广角相机盘左和盘右观测数据分别进行标定解算的结果。由表 3与表 4可知：广角相机单距离观测的像主点坐标存在明显变化，其中u0分量有随摄影距离增加而减小的趋势，v0分量有随摄影距离增加而增大的趋势，体现了偏心元素与相机参数耦合对相机标定的影响；单距离观测相机主距与多距离观测主距最大差异为59 pixel，对应引起图像边缘的角度最大变化约23″；单盘解算平均重投影误差为0.72 pixel，对应图像分辨率15.5″，多距离解算平均重投影误差为0.6 pixel，对应的图像角分辨率约13.3″。
表 3 广角相机单距离观测标定解算Table 3. Calibration solution of wide-angle camera for single distance observationface s/m u0/pixel v0/pixel fx/pixel fy/pixel γ/(″) mr/pixel 0 8.11 1033.24 1267.07 9438.46 9433.63 -64.21 0.91 1 8.11 1030.93 1265.02 9439.43 9434.14 -71.99 0.86 0 12.27 1009.08 1273.66 9414.37 9409.81 -37.40 0.92 1 12.27 1007.02 1271.27 9414.72 9410.60 -46.64 0.88 0 15.83 998.50 1276.43 9407.02 9402.95 -46.28 0.76 1 15.83 996.38 1273.87 9405.88 9401.83 -52.03 0.80 0 19.85 991.33 1278.15 9405.75 9402.92 -55.70 0.57 1 19.85 989.26 1276.18 9407.07 9401.63 -51.05 0.61 0 23.93 986.43 1279.06 9412.55 9411.97 -15.19 0.61 1 23.93 984.45 1277.25 9409.08 9406.26 -49.37 0.47 0 29.70 981.64 1280.41 9396.30 9401.08 -19.98 0.80 1 29.70 979.66 1278.31 9400.80 9401.36 -41.03 0.48 表 4 广角相机多距离观测标定解算Table 4. Calibration solution of wide-angle camera for single distance observationface u0/pixel v0/pixel fx/pixel fy/pixel α/(″) β/(″) γ/(″) x0/mm y0/mm z0/mm k1/10-9 k2/10-15 mr/pixel 0 962.48 1285.71 9380.27 9379.64 -50.16 40.97 -45.15 60.58 -15.88 -48.09 3.33 -3.26 0.61 1 960.61 1283.43 9384.54 9383.22 -50.16 40.97 -50.53 60.22 -15.88 -45.34 2.32 -2.32 0.58 望远镜相机多距离观测数据进行整体观测的重投影误差图如图 5所示。广角相机多距离观测数据进行整体观测的重投影误差图如图 6所示。图 5、图 6为利用相机多距离观测标定结果作为相机参数计算的靶标理论像素坐标与中心提取算法获得的靶标中心像素坐标的误差矢量图。从图中可发现，图像中心视域重投影误差较小，视域边缘重投影误差相对较大，说明相机中心视域内的图像测量精度更高。
![图 5 望远镜相机整体解算重投影误差图]()
图 5 望远镜相机整体解算重投影误差图Figure 5. Overall solution of reprojection error diagram of telescope camera![图 6 广角相机整体解算重投影误差图]()
图 6 广角相机整体解算重投影误差图Figure 6. Overall solution of reprojection error diagram of wide-angle camera从实验结果来看，有：(a)望远镜单距离简化模型中，望远镜相机可视作同轴安置，不同摄影距离相机像主点位置较为稳定，变化量主要体现在分量引起的主距变化；广角相机应视为偏心安置，偏心元素对相机标定解算的影响极大，主要表现为分量引起相机像主点位置在不同摄影距离下明显变化；(b)通过多距离双盘观测联合解算模型可以有效解决标定参数的耦合问题，由此得到的相机参数适用于不同摄影距离的图像测量，其中望远镜相机多距离变焦摄影图像的角度分辨率可达2.0″，单个距离定焦摄影时的角度分辨率可达0.3″；广角相机多距离摄影图像角度分辨率达13.3″；只利用相机的中心视域进行图像测量有更高的角度分辨率，相机标定结果可满足高精度的图像测量要求；(c)望远镜相机多距离多参数模型的标定解算精度低于单距离简化模型标定解算精度，这主要是因为望远镜相机采用了对焦系统，不同距离拍摄过程中需要重新调焦，由此会带来相机焦距和主光轴方向的变化。
在传统相机标定法中，标定场是由全站仪人工测量的系列非共面点的绝对坐标组成，采用这种方式整个控制网的布设过程费时耗力且成本高昂，且控制点的布设还存在通视上的限制；有时囿于加工精度的制约，标定板和标定块不能无限制增大，在广角相机标定时和在较远距离进行相机标定时，标定板和标定块在成像视域内所占图幅较小，控制点难以在成像视域内均匀分布，这会降低传统相机标定解算精度。同时，传统相机标定法要求对控制场进行多距离多视角的摄影，拍摄过程中更多的是使用人工拍摄，拍摄作业的自动化程度低。这些存在的问题，在本文中都得到了克服，很大程度上保证了方法的优势。
4. 结论
采用单目标旋转摄影的图像全站仪相机标定方法，建立了多距离双盘观测的联合解算模型，可以实现广角相机和望远镜相机的精确标定。
(a) 与传统的摄影测量标定方法相比，采用基于旋转摄影的图像全站仪标定方法，无需建立精密的标定控制网，可实现免控制点的精密标定，标定过程可做到自动化。
(b) 多距离双盘观测的联合解算模型能很好地克服参数相关问题，标定解算得到的相机参数可用于不同摄影距离的图像测量，最终望远镜相机图像角度分辨率达2.0″，定焦摄影时图像角度分辨率达0.3″；广角相机图像角度分辨率达13.3″，这个结果优于传统的方法。
(c) 可进一步推广，将相机搭载在全站仪上，采用本文中提出的旋转摄影相机标定方法和解算模型，实现其它相机免控制网的全自动高精度标定。
本文中的研究结果有望提高全站仪的测量精度。
-
![]()
图 4 静水环境系统图
a—工件加工示意图 b—光路示意图
Figure 4. Photo of water environment
a—workpiece processing diagram b—optical path diagram
![]()
图 5 真空环境系统图
a—工件加工示意图 b—光路示意图
Figure 5. Photo of vacuum environment
a—workpiece processing diagram b—optical path diagram
![]()
图 6 空气环境中激光烧蚀截面图
a—模拟结果 b—实验烧蚀表面形貌
Figure 6. Cross section of the laser ablate under air environment
a—simulation results b—experimental results
![]()
图 7 纳秒激光烧蚀单晶硅温度场
Figure 7. Temperature field of nanosecond laser ablation of mono-crystalline silicon
![]()
图 8 纳秒激光烧蚀单晶硅速度场及实验结果
Figure 8. Velocity field of nanosecond laser ablation of mono-crystalline silicon and experimental results
![]()
图 9 不同能量密度烧蚀单晶硅表面形貌
Figure 9. Surface morphology of the ablated mono-crystalline silicon with di-fferent energy densities
-
[1] 王佳敏, 季艳慧, 梁志勇, 等. 532 nm皮秒脉冲激光对单晶硅的损伤特性研究[J]. 中国光学, 2022, 15(2): 242-250. https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGA202202005.htm WANG J M, JI Y H, LIANG Zh Y, et al. Damage characteristics of a 532 nm picosecond pulse laser on monocrystalline silicon[J]. Chin-ese Optics, 2022, 15(2): 242-250(in Chinese). https://www.cnki.com.cn/Article/CJFDTOTAL-ZGGA202202005.htm
[2] 聂世琳, 管迎春. 紫外激光器及其在微加工中的应用[J]. 光电工程, 2017, 44(12): 1169 -1179. https://www.cnki.com.cn/Article/CJFDTOTAL-GDGC201712008.htm NIE Sh L, GUAN Y Ch. Review of UV laser and its applications in micromachining[J]. Opto-Electronic Engineering, 2017, 44(12): 1169-1179(in Chinese). https://www.cnki.com.cn/Article/CJFDTOTAL-GDGC201712008.htm
[3] 卢红. 不同环境下激光与材料的相互作用[D]. 天津: 天津大学, 2017: 10-17. LU H. Interaction between laser and materials under different environmental conditions[D]. Tianjin: Tianjin University, 2017: 10-17 (in Chinese).
[4] 张晋平. 飞秒激光加工高深径比微孔的研究[D]. 上海: 上海交通大学, 2015: 20-26. ZHANG J P. Research of high-aspect ratio holes drilled by femtose-cond laser[D]. Shanghai: Shanghai Jiao Tong University, 2015: 20-26(in Chinese).
[5] UENO I, SHOJI M. Thermal-fluid phenomena induced by nanosecond-pulse heating of materials in water[J]. Journal of Heat Transfer, 2001, 123(6): 1123-1132. DOI: 10.1115/1.1409264
[6] 周翔, 段军, 陈航, 等. 水辅助激光无重铸层钻孔Al2O3陶瓷实验研究[J]. 激光技术, 2018, 42(2): 271-275. DOI: 10.7510/jgjs.issn.1001-3806.2018.02.025 ZHOU X, DUAN J, CHEN H, et al. Experimental study about water-assisted laser drill on Al2O3 ceramics without recast layer[J]. Laser Technology, 2018, 42(2): 271-275(in Chinese). DOI: 10.7510/jgjs.issn.1001-3806.2018.02.025
[7] 朱帅杰, 张朝阳, 储松林, 等. 海量微孔水辅助法皮秒激光加工技术的研究及应用[J]. 中国激光, 2020, 47(3): 03020028. https://www.cnki.com.cn/Article/CJFDTOTAL-JJZZ202003015.htm ZHU Sh J, ZHANG Ch Y, CHU S L, et al. Research and application of massive micropores water-assisted picosecond laser processing technology[J]. Chinese Journal of Lasers, 2020, 47(3): 03020028 (in Chinese). https://www.cnki.com.cn/Article/CJFDTOTAL-JJZZ202003015.htm
[8] DOU H Q, YAO C Z, LIU H, et al. Femtosecond laser ablation of Al-Mg alloy in vacuum and air[J]. Applied Surface Science, 2018, 447(31): 388-392.
[9] WANG Q W, YE Y C, PENG Y, et al. Friction and cutting characteristics of micro-textured diamond tools fabricated with femtosecond laser[J]. Tribology International, 2020, 154: 106720.
[10] 王卫江, 贾凯, 房瑞娜, 等. 激光诱导等离子体辐射特性的研究综述[J]. 激光技术, 2022, 46(4): 499-510. DOI: 10.7510/jgjs.issn.1001-3806.2022.04.010 WANG W J, JIA K, FANG R N, et al. Review on radiation features of laser-induced plasma[J]. Laser Technology, 2022, 46(4): 499-510(in Chinese). DOI: 10.7510/jgjs.issn.1001-3806.2022.04.010
[11] WANG X, QIN Y, WANG B, et al. Numerical and experimental study of the thermal stress of silicon induced by a millisecond laser[J]. Applied Optics, 2011, 50 (21) : 3725-3732. DOI: 10.1364/AO.50.003725
[12] 张廷忠. 毫秒激光打孔过程熔融喷溅、重铸层和微裂纹形成机理研究[D]. 南京: 南京理工大学, 2017: 31-33. ZHANG T Zh. Study on the mechanism of melt ejection recast, layer and micro crack formation in millisecond laser drilling[D]. Nanjing: Nanjing University of Science and Technology, 2017: 31-33(in Chinese).
[13] 钟发成, 吕雪明, 李佳桂, 等. 组合激光对单晶硅热作用的数值分析[J]. 激光技术, 2017, 41(5): 637-643. DOI: 10.7510/jgjs.issn.1001-3806.2017.05.004 ZHONG F Ch, LV X M, LI J G, et al. Numerical analysis of thermal effect of the combined laser on single crystal silicon polymers irradiated by high power laser[J]. Laser Technology, 2017, 41(5): 637-643(in Chinese). DOI: 10.7510/jgjs.issn.1001-3806.2017.05.004
[14] 齐立涛, 李雪, 刘亚升. 单脉冲266 nm固体激光烧蚀蓝宝石的工艺和机理[J]. 黑龙江科技大学学报, 2020, 30(4): 416-421. https://www.cnki.com.cn/Article/CJFDTOTAL-HLJI202004012.htm QI L T, LI X, LIU Y Sh. Process and mechanism for sapphire ablation by singlepulse 266 nm solid-state laser[J]. Journal of Heilongjiang University of Science and Technology, 2020, 30(4): 416-421 (in Chinese). https://www.cnki.com.cn/Article/CJFDTOTAL-HLJI202004012.htm
[15] 曹州权, 谢小柱, 陈蔚芳, 等. 液体辅助激光加工过程中压力检测及应用研究进展[J]. 光电工程, 2017, 44(4): 381-392. https://www.cnki.com.cn/Article/CJFDTOTAL-GDGC201704005.htm CAO Zh Q, XIE X Zh, CHEN W F, et al. Research progress of pressure detection and applications in liquid-assisted laser machining[J]. Opto-Electronic Engineer, 2017, 44(4): 381-392(in Chin-ese). https://www.cnki.com.cn/Article/CJFDTOTAL-GDGC201704005.htm
[16] 齐立涛, 刘凤聪. 266 nm纳秒固体激光在CH薄膜上打孔的工艺实验研究[J]. 激光技术, 2022, 46(6): 767-772. DOI: 10.7510/jgjs.issn.1001-3806.2022.06.009 QI L T, LIU F C. Experimental investigation on 266nm nanosecond solid-state laser drilling on CH film[J]. Laser Technology, 2022, 46(6): 767-772(in Chinese). DOI: 10.7510/jgjs.issn.1001-3806.2022.06.009
[17] 鲍家定. 水层特性对水下激光划切单晶硅质量的影响研究[D]. 桂林: 桂林电子科技大学, 2021: 28-34. BAO J D. Study on the influence of water layer characteristics on quality of under water laser cutting monocrystalline silicon[D]. Gui-lin: Guilin University of Electronic Technology, 2021: 28-34(in Ch-inese).
[18] 彭根琛. 铝合金真空激光焊接特性及能量传输机理研究[D]. 哈尔滨: 哈尔滨工业大学, 2020: 20-25. PENG G Ch. Research on characteristics and energy transfer mechanism of laser welding under vacuum for aluminum alloys[D]. Harbin: Harbin Institute of Technology, 2020: 20-25 (in Chinese).
[19] 张潇允, 张巍, 夏盛强, 等. 高功率激光辐照CFRP的温度场和应力场的数值分析[J]. 激光技术, 2021, 45(5): 636-641. DOI: 10.7510/jgjs.issn.1001-3806.2021.05.017 ZHANG X Y, ZHANG W, XIA Sh Q, et al. Study on numerical analysis of temperature field and stress field of carbon fiber reinforced polymers irradiated by high power laser[J]. Laser Technology, 2021, 45(5): 636-641(in Chinese). DOI: 10.7510/jgjs.issn.1001-3806.2021.05.017
[20] LI G L, SHI D, H H S, et al. Research on the mechanism of laser drilling alumina ceramics in shallow water[J]. The International Journal of Advanced Manufacturing Technology, 2022, 118(11): 3631-3639.
下载:
计量
- 文章访问数: 531
- HTML全文浏览量: 0
- PDF下载量: 7