-
从明暗变化恢复形状,是通过物体的一幅或多幅2D灰度图像获得该物体3D表面形状信息的方法。其原理是利用成像表面每一点亮度值的变化,解析出物体表面的矢量信息,转化为物体表面的深度信息[19]。
对实际图像而言,表面上任意点的亮度受到众多因素影响,如:现实光源、物体表面材料性质和形状,以及相机的位置和参数等。假设物体表面满足朗伯体余弦定律,其表面反射模型如图3所示,物体表面点的图像亮度只取决于该点光源入射角的余弦,即O´点的亮度仅由O点光源入射角的余弦
${\theta _i}$ 决定[16]。为了进行表面形貌的重建,需要对拍摄图像中的灰度信息进行转化。基于Alhulayil等人的研究[20],此时表面反射模型可变换为:
$$ \begin{split} E\left( {x,y} \right) =& R\left( {p\left( {x,y} \right),q\left( {x,y} \right)} \right) =\\ & {I_0}\rho \dfrac{{{p_0}p + {q_0}q + 1}}{{\sqrt {p_0^2 + q_0^2 + 1} \sqrt {p_{}^2 + q_{}^2 + 1} }} \end{split} $$ (1) 式中:
$E\left( {x,y} \right)$ 为图像亮度;$R\left( {p,q} \right)$ 为反射函数;I0为入射光源强度;ρ为物体表面反射率;p、q为表面高度沿x、y方向的变化率;p0、q0由入射光源方向决定。如公式(1)所示,此时物体表面存在p、q两个未知量,而图像灰度信息仅能提供一个约束。一般情况下,由该模型所确定的SFS问题是病态的,无法求出唯一解[21]。以上模型尚无法满足烧蚀坑表面形貌重建的重构要求。因此,需要添加约束条件对重建模型进行约束。为了进行烧蚀坑表面形貌的重建,需要物体表面的信息对图像灰度进行分析。为了构建物体目标表面与图像灰度表面之间的联系,采用了基于灰度梯度约束的思路进行烧蚀坑的表面形貌重建[16],其原理如图4所示。
假设光源为垂直入射,此时光源的向量为
$L( {0,} $ $ {0,1} )$ ,以镜头中心为原点坐标。坐标系的z轴与成像系统的光轴重合,并且x轴和y轴平行于成像平面。在成像过程中,将物体表面投影到图像平面时,将表面高度转化为灰度值。此时,Ai为图像平面上的任意点;Ag是法向量为Ng(−pg,−qg,1)且倾斜角度为θg时,其在灰度平面上的对应点;而Aw是法向量为N(−p,−q,1)且倾斜角度为θw时,其在目标表面上的对应点。
$\varphi $ 为以上点所对应的偏角。为了约束灰度表面方向与目标表面方向梯度的关系,由解析几何理论可知[16],Ag与Aw在图像平面上对应的偏角均为
$\varphi $ 。因此,目标表面和灰度表面间的梯度关系可表示为:$$\frac{q}{p} = \frac{{{q_g}}}{{{p_g}}} = \tan \varphi $$ (2) 此时物体表面梯度定义为:
$$\left\{ {\begin{array}{*{20}{c}} {p = k{p_g}} \\ {q = k{q_g}} \end{array}} \right.$$ (3) 式中:k为梯度比例系数。将其作为表面形貌重建问题的约束条件,此时灰度梯度(pg,qg)可由图像
$E\left( {x,y} \right)$ 的灰度给出:$$\left\{ {\begin{array}{*{20}{c}} {{p_g} = E{{(x,y)}_x}} \\ {{q_g} = E{{(x,y)}_y}} \end{array}} \right.$$ (4) 其中,
$E\left( {x,y} \right)$ x和$E\left( {x,y} \right)$ y分别表示沿x、y方向的图像灰度梯度。因此物体表面梯度可写为:$$ \left\{\begin{array}{c}p=k \cdot E{(x,y)}_{x}\\ q=k \cdot E{(x,y)}_{y}\end{array}\right.$$ (5) 此时任意位置的图像灰度梯度可由相邻两点的灰度值的差值进行计算:
$$\left\{ {\begin{array}{*{20}{c}} {E{{(x,y)}_x} = (E(x + 1,y) - E(x - 1,y))/2} \\ {E{{(x,y)}_y} = (E(x,y + 1) - E(x,y - 1))/2} \end{array}} \right.$$ (6) 由入射光源的向量
$L\left( {0,0,1} \right)$ 可知,向量p0、q0均为0,将公式(5)代入公式(1)时,可得:$$E(x,y) = \frac{{{I_0}\rho }}{{\sqrt {1 + {k^2}\left[ {E(x,y)_x^2 + E(x,y)_y^2} \right]} }}$$ (7) 此时梯度比例系数为:
$$k(x,y) = \pm {\left( {\frac{{{{\left( {{I_0}\rho } \right)}^2} - E{{(x,y)}^2}}}{{E{{(x,y)}^2}\left( {E(x,y)_x^2 + E(x,y)_y^2} \right)}}} \right)^{1/2}}$$ (8) 将公式(8)代入公式(5)后,物体表面梯度可表示为:
$$ \left\{\begin{array}{c}p=\pm {\left(\dfrac{{\left({I}_{0}\rho \right)}^{2}-E{(x,y)}^{2}}{E{(x,y)}^{2}\left(E{(x,y)}_{x}^{2}+E{(x,y)}_{y}^{2}\right)}\right)}^{1/2} \cdot E{(x,y)}_{x}\\ q=\pm {\left(\dfrac{{\left({I}_{0}\rho \right)}^{2}-E{(x,y)}^{2}}{E{(x,y)}^{2}\left(E{(x,y)}_{x}^{2}+E{(x,y)}_{y}^{2}\right)}\right)}^{1/2} \cdot E{(x,y)}_{y}\end{array}\right.$$ (9) 其中,正负号仅能影响重建形状的凸凹,对相对高度的具体数值计算无影响。此时可通过公式(9)得到物体表面任意位置的表面梯度p、q,从而进行物体表面梯度信息和物体表面高度信息的转化。
假设z0为图像中初始点
$\left( {{x_0},{y_0}} \right)$ 的高度,物体表面任意点相对高度的计算公式可表示为:$${\textit{z}}\left( {x,y} \right) = {{\textit{z}}_0} + \int_{\left( {{x_0},{y_0}} \right)}^{\left( {x,y} \right)} {p{\rm{d}}x + q\rm{d}}y $$ (10) 此时通过公式(10)可进行物体表面高度的计算。
-
根据上述重构原理,重构表面形貌时首先要确定参数I0和参数ρ。而在SEM的真空测量环境中,两参数的具体数值无法直接得到。因此,文中基于朗伯体表面的反射特性对计算过程进行简化:从任意角度观察到的物体表面上同一点的亮度相同,即物体表面各点的反射光强与光源入射角的余弦成正比[16, 20-21]。此时公式(1)可转化为:
$$E\left( {x,y} \right) = {I_0}\rho \cos {\theta _i}$$ (11) 式中:
${\theta _i}$ 为入射光与表面法向量间的夹角;I0ρ为入射光源强度和物体表面反射率相乘后得到的比例常数。在入射光源强度一定的情况下,由于工质侧面任意位置处表面反射率相等,此时二者相乘后得到的常数参数为恒定值。为了进行烧蚀坑表面形貌的重建,文中给出了一种参数计算方法。
(1)参数I0ρ的计算
在SFS重建模型的应用中[16, 20-21],入射光源强度I0和物体表面反射率ρ可分别由仪器直接测量和查找反射率表的方式得到,不涉及参数I0ρ的求解问题。而在该研究中,对百微米级尺度烧蚀坑的侧面图像进行拍摄时,使用了扫描电子显微镜(SEM)进行拍摄。由于SEM需要在真空环境中进行图像拍摄,此时入射光源强度I0无法直接测量;而对于研究过程中使用的GAP含能工质,其物体表面反射率ρ同样无法直接得到。因此,参数I0ρ的求解成为了SFS方法应用时所关注的问题。
在单脉冲激光烧蚀情况下,假设的烧蚀坑形状和真实烧蚀坑形状通常如图5所示。其中,图5(a)为假设情况下的烧蚀坑形状,Z0、Z1为烧蚀坑侧面图像上已知相对高度和距离的两点,S1、S2分别为烧蚀坑上下表面轮廓,其理想烧蚀形状如图5(b)所示;而在实际烧蚀中,烧蚀坑侧面通常如图5(c)所示,由于各种随机因素的影响,常为不规则形状。
图 5 烧蚀坑形状。(a)假设烧蚀坑形状; (b)真实烧蚀坑上下表面; (c)真实烧蚀坑侧面
Figure 5. Ablation pit shape. (a)Assuming the shape of the ablation pit;(b) Upper and lower surfaces of real ablation pits; (c) The side of the real ablation pit
根据图5(b)中SEM得到的烧蚀坑图像,可以得到烧蚀坑上下表面的具体大小和尺寸,在上下表面的同一方向上,任意两点Z0、Z1的相对高度h均可通过图像测量得到。
因此,由公式(10)可知,两点间相对高度为:
$${{\textit{z}}_1} = {{\textit{z}}_0} + \int_{{y_0}}^{{y_1}} {q{\rm{d}}y} $$ (12) 在实际情况中,需要通过物体表面梯度值和初始点的假定高度对某点的相对高度进行计算。将公式(9)代入公式(11)中,可化为:
$${{\textit{z}}_1} = {{\textit{z}}_0} + \int_{{y_0}}^{{y_1}} { \pm {{\left\{ {\dfrac{{{{\left( {{I_0}\rho } \right)}^2} - E{{\left( {x,y} \right)}^2}}}{{E{{\left( {x,y} \right)}^2}\left[ {E\left( {x,y} \right)_x^2 + E\left( {x,y} \right)_y^2} \right]}}} \right\}}^{\frac{1}{2}}}E{{\left( {x,y} \right)}_y}{\rm{d}}y} $$ (13) 其中,对于同一幅图像,未知量参数I0ρ可视为常数,其余参数均可由图像信息得到。公式(14)中对参数I0ρ的求解过程过于繁杂,且难以进行计算。因此,文中通过采用了迭代的思想对参数进行计算,此时函数
$h\left( {{I_0}\rho } \right)$ 如公式(14)所示。计算中,在设置参数的初始值后,不断进行迭代,当
$h\left( {{I_0}\rho } \right)$ 的值小于0.001时,认为得到参数I0ρ的数值。此时可通过参数I0ρ进行烧蚀坑表面形貌的重建。$$\begin{array}{l} h\left( {{I_0}\rho } \right) = \left( {{{\textit{z}}_1} - {{\textit{z}}_0}} \right) -\\ \begin{array}{*{20}{c}} {}&{} \end{array} \displaystyle\int_{{y_0}}^{{y_1}} { \pm {{\left\{ {\frac{{{{\left( {{I_0}\rho } \right)}^2} - E{{\left( {x,y} \right)}^2}}}{{E{{\left( {x,y} \right)}^2}\left[ {E\left( {x,y} \right)_x^2 + E\left( {x,y} \right)_y^2} \right]}}} \right\}}^{\frac{1}{2}}}E{{\left( {x,y} \right)}_y}{\rm{d}}y} \\ \end{array} $$ (14) (2)噪声的平滑处理
在烧蚀坑表面形貌重建过程中,灰度值图像通常是由SEM得到的烧蚀坑侧面图像对连续变化的灰度值进行采样和量化获得的[22]。图像中会出现不均匀的灰度噪声,导致灰度值图像中烧蚀坑内部区域存在随机误差。由于在图像水平方向上,以烧蚀坑边缘两侧为边界时,发现曲线灰度值存在明显的差异。采用分段拟合的方式,通过边缘检测得到边缘点,对这些点进行曲线拟合后得到烧蚀坑两侧的边界轮廓曲线,此时将拟合后得到的曲线作为烧蚀坑图像分段的依据,通过三次多项式拟合函数对各段灰度值进行曲线拟合,从而对各分段烧蚀坑灰度值进行了有效地平滑。
如图5(b)所示,由于在图像水平方向上的灰度值存在明显的差异,在进行灰度图像分段,拟合灰度值曲线时可能会出现跳跃间断点。为保证重建曲线的平滑和重建表面形貌的准确,同时尽可能减小重建边界处的不稳定性,需要结合分段点两侧灰度值对分段点进行赋值,从而对分段点进行约束。此时,通过约束分段点可得到平滑后图像灰度值曲线,进而对灰度值图像进行表面形貌重建。
-
基于以上分析,烧蚀质量计算的算法原理如图6所示。该算法主要包括5个步骤:原始图像的放缩、灰度值图像的平滑和降噪、约束条件下表面梯度值的计算、物体表面相对高度计算以及烧蚀坑体积和总烧蚀质量的估算。
(1)读取原始图像数据。结合真实尺寸与像素数量的关系,对图像进行放缩,使图像中任意两点间的真实距离和像素间距为1∶1;
(2)得到放缩后图像的灰度值图像
${E_0}(x,y)$ 后,为了平滑降噪,通过提取边缘轮廓曲线作为各方向上的分段依据,进行分段曲线拟合。考虑到分段点处可能存在跳跃间断点,需要结合分段点两侧灰度值对分段点进行赋值,以此对分段点处灰度值进行约束,最终得到处理后灰度值图像$E(x,y)$ ;(3)基于处理后的灰度值图像
$E(x,y)$ ,由公式(6)可得到相应的灰度梯度$({E_x},{E_y})$ 。此时可通过公式(9)计算图像中每点灰度值$E(x,y)$ 对应的物体表面梯度p、q,式中参数I0ρ可通过图5(b)、(c)中已知相对位置和高度的两点,通过迭代的思想,采用公式(14)进行计算。(4)将已知具体位置和高度的点
$\left( {{x_0},{y_0}} \right)$ 设置为初始点,其初始高度为z0。因此,物体表面上任意一点的相对高度可采用公式(10)进行计算,此时可得到重建后烧蚀坑任意位置处深度及烧蚀坑的表面形貌。(5)根据相对高度进行烧蚀坑体积的估算。对剖开的两侧烧蚀坑分别利用算法进行体积的估算,并进行求和。此时单脉冲烧蚀质量计算公式为:
$$M = {\rho _0}\left( {{V_1} + {V_2}} \right)$$ (15) 式中:M为单脉冲烧蚀质量;ρ0为烧蚀工质密度;V1和V2分别为两侧重建烧蚀坑体积。
-
为了对重建算法的有效性和准确性进行验证,研究中对合成图像和烧蚀坑图像进行了实验。
-
基于算法研究中的基准图像[23-24],以及激光烧蚀的模拟烧蚀坑形状[6],实验中采用了半圆台的图像,由曲面合成的图像如图7(a)所示,上表面直径为400 pixel,下表面直径为900 pixel,高度为1 000 pixel。此时假设合成图像符合朗伯体反射模型,其中光源矢量方向为(0,0,1)。
图 7 半圆台曲面的仿真结果。(a) 半圆台曲面 ;(b)合成图像 ;(c)重建曲面 ;(d)重建误差
Figure 7. Simulation results of semi-circular surface. (a) Semicircle surface; (b) Composite image; (c) Reconstruction surface; (d) Reconstruction error
图7(b)为原始曲面的合成图像,其重建表面见图7(c),图7(d)显示了原始曲面和重构曲面间的误差。可看出重构曲面与原始曲面相似,最大误差值为137 pixel。重建误差主要发生在圆台边界附近,其他部分误差可忽略不计。这是由于边界处灰度梯度变化的不连续性,使得边界附近梯度变化剧烈,最终导致了重建灰度值的异常,使得重建结果出现误差。
为了评价重建算法的性能,利用平均绝对误差(MAE)和均方根误差(RMSE)对原始曲面和重建曲面进行了比较,结果如表1所示。
表 1 MAE和RMSE的比较结果
Table 1. Comparison results of MAE and RMSE
Composograph MAE/pixel RMSE/pixel Circular truncated cone 30 45 表1中,重建表面与原始曲面间的MAE为30 pixel、最大RMSE仅为45 pixel。结果表明,重构后的曲面与原始曲面具有较好的一致性,重建算法具有良好的重建效果和较高的实验精度。
-
在实验中通过SEM收集,以工作电流15 A,脉宽400 μs状态下的半导体激光器,烧蚀基底为50 μm厚的PI材料,工质为200 μm厚,密度为1 600 kg/m3的掺杂3%纳米碳粉GAP工质时,尺寸为1 024×943 pixel的一侧单脉冲烧蚀坑侧面图像。
由于物体表面真实尺寸和图像尺寸间存在差异,为了解决重建物体表面形貌时的尺寸问题,在对图像使用重建算法前进行放缩,使两者尺寸大小为1∶1。此时处理后的典型侧面烧蚀坑图像如图8(a)所示,通过公式(14)迭代计算得到的中间参数I0ρ值为155。
图 8 烧蚀坑图像。(a)侧面图;(b)传统体积法;(c)重建体积法;(d)重建图;(e)重建俯视图
Figure 8. Image of ablation pit. (a) A profile view; (b) Traditional volume method; (c) Reconstruction volume method; (d) Reconstruction plan; (e) Reconstruction of the top view
当激光与工质相互作用时,靠近透明基底侧区域较早开始吸收激光能量,形成烧蚀区。由于GAP工质掺杂了纳米碳粉,碳粉在工质内形成团聚区域,在吸收激光能量形成局部高温后,导致其周围的聚合物被迅速气化和离化[25]。随着冲击波的形成以及激光烧蚀推进过程的进行,烧蚀区压力得到释放,区域温度逐渐降低,烧蚀坑内部及其周围的工质以凹凸不平的状态凝结。此时,在侧面图像中可以发现明显的不规则区域。在获取图像信息时,这些区域会导致信息存在一定程度的误差。为了削弱不规则区域对重建造成的影响,保证重建形貌与烧蚀坑轮廓吻合,在进行烧蚀坑表面形貌重建时,需要对图像进行平滑处理。
如图8(b)所示,传统体积估算方法中,通常在已知烧蚀坑上下表面形状的情况下,采用圆柱或圆台形状对烧蚀坑形貌进行模拟[6]。然而从图像中可以发现,真实烧蚀坑形状与模拟形状之间存在差异。即真实烧蚀坑通常为不规则形状,其体积小于模拟形状体积,且对应烧蚀质量小于模拟烧蚀质量。由于计算模型与真实情况存在差异,会导致计算的推进性能参数与真实值存在偏差。为此,文中采用了如图8(c)所示的边缘轮廓作为烧蚀坑边界,通过表面形貌重建算法对不规则烧蚀坑进行体积的重构,从而进行烧蚀质量的计算。采用传统体积法与表面形貌重建方法计算的体积和质量如表2所示,后者计算的数值小于前者,且更接近真实烧蚀坑的数值。相较于计算模型无法准确反映烧蚀坑形状的传统体积方法,表面形貌重建方法计算的烧蚀坑质量更为接近真实结果。
表 2 重建烧蚀坑参数
Table 2. Reconstructed ablation pit parameters
Ablation pit parameters Volume/μm3 Quality/μg Ideal frustum of a cone 7.54×106 12.07 Reconstruction of ablation pits 4.45×106 7.12 由于成像过程中不可避免地受到二次反射、噪声或杂散辐射的干扰,使得成像结果在某些细节上存在差异。此时经算法重构后,不同方向上的物体表面形貌分别如图8(d)、(e)所示。从重建烧蚀坑的大致形貌可发现,用该方法重建的烧蚀坑已基本满足烧蚀坑形状重建的需求。
为了对重建方法的可重复性进行检测,采用如图9(a)所示的另一幅烧蚀坑图像进行烧蚀坑表面重建。图9为工作电流25 A,脉宽400 µs状态下的半导体激光器,通过烧蚀同一靶带的不同位置得到的,其重建结果如图9(d)、(e)所示,通过公式(14)迭代计算得到的中间参数I0ρ值为100,其对照烧蚀坑体积和质量如表3所示。
图 9 对照烧蚀坑图像。(a)侧面图;(b)传统体积法;(c)重建体积法;(d)重建图;(e)重建俯视图
Figure 9. Contrast ablation pit image. (a) A profile view; (b) Traditional volume method; (c) Reconstruction volume method; (d) Reconstruction plan; (e) Reconstruction of the top view
表 3 对照烧蚀坑参数
Table 3. Contrast ablation pit parameters
Ablation pit parameters Volume/µm3 Quality/µg Ideal frustum of a cone 6.27×106 10.03 Reconstruction of ablation pits 5.34×106 8.544 如图8(d)和图9(d)所示,此时利用SFS方法重建的烧蚀坑与图8(a)和图9(a)所示的原始烧蚀坑相比,前后形状基本吻合,得到了较为理想的表面形貌恢复效果。而如表2和表3所示,尽管脉冲激光烧蚀工质时,烧蚀坑形状存在差异,但通过重建算法的两次单脉冲激光烧蚀坑进行计算,可以发现两次得到的烧蚀坑体积与质量,均与假设其为理想圆台形状时的计算结果处于同一数量级,且小于圆台形状的计算数值。当以相同条件重复烧蚀100次时,总烧蚀质量可达到mg量级,与罗乐乐等人的研究结果相符[26]。因此,在对重建烧蚀坑的体积与质量进行计算时,可以得到较为合理精确的计算结果。此时可以通过该方法重建后的烧蚀坑图像,进行百微米级烧蚀体积和质量的估算。
Estimation method of laser micro-scale ablation quality based on image
-
摘要: 激光微尺度烧蚀下的烧蚀质量是激光微推进系统中评估比冲等性能参数的重要物理量,文中在利用扫描电镜得到的烧蚀坑图像的基础上,基于灰度表面重构形貌方法,建立了微尺度烧蚀坑形貌三维模型重构算法;针对图像噪声影响,对图像进行分段平滑,有效地解决了量化误差和灰度模糊对重构形貌的影响,实现了烧蚀坑体积估算。结果表明:所提出的方法能够快速有效地重构激光微尺度烧蚀下烧蚀坑形貌的三维形状,进而计算其烧蚀质量。Abstract: The ablation quality under laser micro-scale ablation is an important physical quantity for evaluating performance parameters such as specific impulse in the laser micro-propulsion system. Based on the ablation pit image obtained by the scanning electron microscope, a micro-scale ablation pit shape three-dimensional model reconstruction algorithm was established based on the gray-scale surface reconstruction method; In view of the effect of image noise, the image was segmented and smoothed, which effectively solved the effect of quantization error and grayscale blur on the reconstructed shape, and realizeing the ablation pit volume estimation. The results show that the proposed method can quickly and effectively reconstruct the three-dimensional shape of ablation pits under laser micro-scale ablation, and then calculate the ablation quality.
-
Key words:
- laser micro-propulsion /
- laser micro-ablation /
- ablation quality /
- SEM image /
- SFS method
-
表 1 MAE和RMSE的比较结果
Table 1. Comparison results of MAE and RMSE
Composograph MAE/pixel RMSE/pixel Circular truncated cone 30 45 表 2 重建烧蚀坑参数
Table 2. Reconstructed ablation pit parameters
Ablation pit parameters Volume/μm3 Quality/μg Ideal frustum of a cone 7.54×106 12.07 Reconstruction of ablation pits 4.45×106 7.12 表 3 对照烧蚀坑参数
Table 3. Contrast ablation pit parameters
Ablation pit parameters Volume/µm3 Quality/µg Ideal frustum of a cone 6.27×106 10.03 Reconstruction of ablation pits 5.34×106 8.544 -
[1] Phipps C R, Luke J R, Helgeson W, et al. Performance test results for the laser powered microthruster[C]//American Institute of Physics, 2006. [2] Shen S, Jin X, Hao C. Cleaning space debris with a space-based laser system [J]. Chinese Journal of Aeronautics, 2014, 27(4): 805-811. doi: 10.1016/j.cja.2014.05.002 [3] Phipps C R. A laser-optical system to re-enter or lower low Earth orbit space debris [J]. Acta Astronautica, 2014, 93: 418-429. doi: 10.1016/j.actaastro.2013.07.031 [4] Phipps C R, Baker K L, Libby S B, et al. Removing orbital debris with lasers [J]. Advances in Space Research, 2012, 49(9): 1283-1300. doi: 10.1016/j.asr.2012.02.003 [5] Esmiller B, Jacquelard C, Eckel H A, et al. Space debris removal by ground-based lasers: main conclusions of the European project CLEANS PACE [J]. Applied Optics, 2014, 53(31): I45-I54. doi: 10.1364/AO.53.000I45 [6] Ye Jifei, Hong Yanji, Wang Guangyu, et al. Measurement method for microgram amounts laser ablation mass [J]. Journal of Propulsion Technology, 2010, 31(3): 372-376. (in Chinese) [7] Lilly T C, Ketsdever A, Pancotti A P, et al. Development of a specific impulse balance for capillary discharge pulsed plasma thrusters [J]. Journal of Pro-pulsion and Power, 2009, 25(3): 823-826. doi: 10.2514/1.40261 [8] Jin Xing, Wu Wentang, Zhou Weijing, et al. Design method of micro-impulse measuring system based on pivots for pulsed micro-thruster [J]. Infrared and Laser Engineering, 2019, 48(S1): S117005. (in Chinese) [9] Jin Xing; Chang Hao; Ye Jifei. Methods of extreme short pulse laser ablation impulse coupling measurem ent [J]. Infrared and Laser Engineering, 2017, 46(3): 0329002. (in Chinese) [10] Chang Hao, Ye Jifei, Lin Zhengguo. Simulation analysis of response characteristics of torsion pendulum by pulse micro-thrust [J]. Infrared and Laser Engineering, 2016, 45(S2): S218001. (in Chinese) [11] Ketsdever A D, D'Souza B C, Lee R H. Thrust stand micromass balance for the direct measurement of specific impulse [J]. Journal of Propulsion and Power, 2008, 24(6): 1376-1381. doi: 10.2514/1.36921 [12] Zhou Weijing, Hong Yanji, Ye Jifei, et al. Direct measurement method of specific im-pulse for pulse laser ablation micropropulsion [J]. Journal of Propulsion Technology, 2017, 38(6): 1434-1440. (in Chinese) [13] Gonzales D A, Baker R P. Micropropulsion using a Nd: YAG microchip laser[C]//SPIE, 2002, 4760: 752-765. [14] Heise B, Stifter D. Quantitative phase reconstruction for ortho gonal scanning differential phase contrast optical coh-erence tomo graphy [J]. Optics Letters, 2009, 34(9): 1306-1308. doi: 10.1364/OL.34.001306 [15] Hu J-F, Zheng W-S, Xie X, et al. Sparse transfer for facial shape- from-shading [J]. Pattern Recognition, 2017, 68: 272-285. doi: 10.1016/j.patcog.2017.03.029 [16] Shi W C, Zheng J M, Li Y, et al. Three-dimensional reconstruction method for machined surface topography based on gray gradient constraints [J]. Applied Sciences, 2019, 9(3): 591. doi: 10.3390/app9030591 [17] Ye Jifei, Hong Yanji, Wang Guangyu, et al. Research progre-ss in micro-laser plasma propulsion [J]. Chinese Optics, 2011, 4(4): 319-326. (in Chinese) doi: 10.3969/j.issn.2095-1531.2011.04.003 [18] Zhang Yi, Zhang Zhiyuan, Lu Xin, et al. Development of laser plasma micro-propulsion [J]. Science & Technology Review, 2007(19): 70-74. (in Chinese) doi: 10.3321/j.issn:1000-7857.2007.19.014 [19] Tozza S, Falcone M. Analysis and approximation of some shape -from-shading models for non- Lam bertian su-rfaces [J]. Journal of Mathematical Imaging and Vision, 2016, 55(2): 153-178. doi: 10.1007/s10851-016-0636-x [20] Alhulayil M, Younes A B, Turner J D. Higher order algorithm for solving La-mbert’s problem [J]. The Journal of the Astronautical Sciences, 2018, 65(4): 400-422. doi: 10.1007/s40295-018-0137-9 [21] Wang G, Cheng J. Three-dime nsional reconstruction of hybrid surfaces using perspective shape from shading [J]. Optik, 2016, 127(19): 7740-7751. [22] Meng Q, Wang H, Xu W, et al. A coupling method incorporating digital image processing and discrete element method for modeling of geomaterials [J]. Engineering Computations, 2018, 35(1): 411-431. [23] Yang L, Zhang N, Ren Y, et al. Linear perspective shape-from-shading method with two images [J]. Journal of Systems Engineering and Electronics, 2015, 26(5): 1080-1087. doi: 10.1109/JSEE.2015.00117 [24] Bruvoll S, Reimers M. Spherical surface parameterization for perspective shape from shading [J]. Pattern Recognition Letters, 2012, 33(1): 33-40. doi: 10.1016/j.patrec.2011.09.026 [25] Lukas Urech, Marc Hauer, Thomas Lippert, et al. Designed polymers for laser- based microthrusters: correlation of thrust with material, plasma, and shockwave properties[C]// SPIE, 2005, 5448: 52-64. [26] Luo Lele, Dou Zhiguo, Zhou Weijin, et al. Effects of different dopants on laser ablation propulsion performance of GAP [J]. Journal of Propulsion Technology, 2018, 39(11): 2615-2623. (in Chinese)