1.   引　言

近年来，激光技术发展迅速，由于激光器具有定向发光、亮度高、能量密度大，相干性和单色性好的特点，已经被广泛应用在医学、工业和军事领域^[1]。激光的功率密度分布、能量分布等关系到激光系统的光束质量、到靶能力、瞄准精度、质心漂移等重要参数，是分析激光大气传输效应和评价激光系统出光性能的技术依据^[2-4]。而在实际应用中，由于受到大气环境和天气情况等因素的影响,激光在传播过程中，其能量密度分布往往不均匀,甚至会出现多个光强峰值。因此，在工业加工、激光传输及目标损伤效应试验等激光应用领域均需要对激光能量分布等参数进行准确测量。

目前，常用的激光强度分布测量方法主要有散射成像法^[5]、扫描取样法^[6]和阵列探测法^[7]等，这些方法具有各自的应用特点，在激光性能评估中发挥着重要作用。然而，上述方法也存在一些缺陷，可以总结为以下三点：第一，在测量时需要对激光束进行空间分束取样，这或多或少会对激光的分束分布参数产生影响，难以做到绝对的非接触式测量；第二，这些方法均为光束空间传输某横截面的光强二维分布测量，无法在单次测量中获得光束传输方向纵向截面的光强分布，进而得到光束焦斑的光强分布及位置信息；第三，测量设备复杂，测试精度不高。

与上述测量方法相比，三维成像技术具有非接触性，精度高，能直接实现三维测量的优势。通常，三维成像的方法主要有扫描平面成像^[8-9]和计算成像。而计算成像又可以分为光场成像^[10-12]和计算层析成像^[13-15]。扫描平面成像通过传递激光片，逐层地扫描被测目标，并使用相机对每一个镜头的信号进行连续的采集。但由于来自不同扫描层的信号是在不同时刻采集的，因此该方法最主要的缺点是不能表示瞬时的三维光场分布。光场成像方法只需一台高速相机即可实现瞬时的三维成像，但是该方法的实现需要权衡空间分辨率和视场数量^[16]，且对于CCD芯片的分辨率要求较高。而计算层析成像技术可以克服上述两种方法的不足，即通过使用多个相机同时采集激发光或者自发光的投影^[17-18]。

本文提出一种基于层析成像技术的激光强度分布测量方法，通过数值模拟验证了该方法的可行性，并加以实验验证。模拟研究中，通过量化不同激光光强分布模型以及不同信噪比下的重建精度，对激光强度分布测量技术的可行性进行评估，即验证成像模型的正确性、重建算法的收敛性及在不同信噪比下的鲁棒性。搭建层析系统进行实验验证，通过测量激光穿过若丹明-乙醇溶液后的体激光诱导荧光信号，并对其进行经过吸收矫正的三维重建，从而间接得到激光光强三维分布。本文提出的方法具有层析成像技术非接触、精度高、可实现三维测量的特点，同时借由内窥层析技术解决了传统层析成像设备造价高昂的问题，在激光应用领域具有广泛前景。

2.   层析成像原理

2.1   层析成像数学建模

如图1所示，感兴趣区域（VOI）在空间中的能量分布可以看作是其中所有点光源的集合。因此，任意一个相机在拍摄时所得到的图像是所有点光源经成像系统在CCD上所成的像的叠加。单位辐射强度的点光源所成的像即为点扩散函数。因此，光强在CCD上成的像即为点扩散函数的权重和，其中权重为空间中各个点光源的强度^[19]。当不考虑自吸收和散射时，某个像素$\left(x_{{\rm{p}}}, y_{{\rm{p}}}\right)$的光强${{p}}({x_{\rm{p}}},{y_{\rm{p}}})$在数学上可以表示为：

图  1  成像过程原理图

Figure  1.  Illustration of the imaging process

下载: 全尺寸图片 幻灯片

其中：${{I}}(x,y,{\textit{z}})$是空间中某一点$\left(x_{i}, y_{i}, {\textit{z}}_{i}\right)$的光强；PSF为点扩散函数。请注意$\left(x_{{\rm{p}}}, y_{{\rm{p}}}\right)$和$\left(x_{i}, y_{i}, {\textit{z}}_{i}\right)$分别是图像坐标系和世界坐标系中的点。如果将空间中的重建区域划分为$N_{\text {voxel }}=\left(N_{x} \times N_{y} \times N_{x}\right)$个体素，则公式（1）可以离散化为：

其中：$I_{i}$表示第i个体素的光强；${{PS}}{{{F}}_{{i}}}({x_p},{y_p})$为第i个体素在CCD上所对应的点扩散函数在某个像素中心的值；$\Delta x$,$\Delta y$,$\Delta \textit{z}$分别为沿着3个坐标轴方向体素的大小。从公式（2）可以看出，实际上像平面上的每一个像素满足一个线性方程，该方程的变量为空间中所有体素的发光光强。因此，对于每一个像素都能得到类似公式（2）的一个方程，该CCD上的所有像素均可构成一个方程组。当采用多个相机在不同角度拍摄时则可获得多个线性方程组。把这些方程组放到一起则可以获得一个总的线性方程组：

其中：$\vec{{{p}}}$为向量，表示所有像素所获得的信号集合；矩阵${{{\boldsymbol{A}}}}$中每一列表示某个体素在CCD上的点扩散函数；$\vec{{{f}}}$则为空间中所有体素的光强。

由公式（3）可知，如要求得$\vec{{{f}}}$则需知道空间中各个体素在CCD上的点扩散函数。首先，需要通过标定的方法测得所有体素和各个相机之间的相对位置^{[20, 21]}。之后，采用蒙特卡洛模拟的方法计算点扩散函数：即在每一个体素中随机产生大量光子，随后，追踪每个光子经过成像系统后落在CCD上的像素位置^[22],记录这个像素的位置。这些光子在CCD上的成像会有一个分布。那么，对于某一个特定的像素而言，可以得到相关体素的贡献关系，也就可以得到点扩散函数。之后便可通过求解公式（3）获得待测三维空间中的待测强度分布。

在对比多种重建算法之后，本工作采用代数重建算法^[23]对上述的反演问题进行求解。该算法在投影数量较少时的表现尤为突出，同时能够有效抑制伪影问题，其迭代方式数学表达为：

其中，$\vec{{{f}}}$代表每一次迭代的方程解，上标n代表迭代次数。${{{{\boldsymbol{A}}}}_{{i}}}$代表系数矩阵${{{\boldsymbol{A}}}}$的第i行，p_i代表投影$\vec{{{p}}}_{i}$的第i个元素，矩阵上标T代表矩阵的转置。λ_ART是松弛因子，用于控制迭代过程的收敛速度，本工作中取值为0.07。

本工作所有有关重建过程的正向建模以及反向求解均基于Visual Studio平台，采取Fortran编程语言。

3.   数值模拟研究

为验证本课题所提出的激光能量分布测量技术的可行性（成像模型的正确性以及重建算法的收敛性及在不同信噪比下的鲁棒性），本工作从以下两方面进行了数值仿真。

3.1   激光光强三维分布对重建精度的影响

首先，为了验证成像模型及重建算法在不同形式激光光强三维分布下的适用性，人为设计了不同的激光光强三维分布模型（phantom）用于数值仿真，包括：激光光强三维分布呈高斯分布以及非高斯分布的phantom（假设沿着激光传播方向任一截面上时激光光强分布不变），以及在实验过程中可能会出现的光束沿着传播方向发散的情况（在每一截面上仍为高斯分布）。传感器均布于水平面内的半圆周上，等角间距20°，距离重建体中心160 mm。基于公式（3）建立重建问题的正向模型，即计算模拟投影：根据采取的三维光强模型以及计算出的点扩散函数，对重建区域内每一体素内的光强求权重和（权重即为点扩散函数），即可得到每一个角度的模拟投影。并基于公式（4）采取代数重建算法进行反向问题求解，即三维重建。在重建区域离散程度均相同的情况下，激光光强的三维分布phantom及重建结果的等值面示意图如图2（彩图见期刊电子版）所示，其中，图2（a）~图2（b）为沿着光束传播方向光强无变化的高斯分布，图2（c）~图2（d）为光束沿传播方向发散的分布（但每一截面上仍为高斯分布），图2（e）~图2（f）为沿着光束传播方向光强无变化的非高斯分布；其中：（a）、（c）和（e）为原始phantom,（b）、（d）和（f）为对应的重建结果。

图  2  不同激光光强三维分布phantom及其重建结果等值面示意图：（a−b）为沿着光束传播方向光强无变化的高斯分布；（c−d）为沿着光束传播方向光强分布有变化的分布；（e−f）为沿着光束传播方向光强无变化的非高斯分布；（a）、（c）、（e）为原始phantom,（b）、（d）、（f）为对应的重建结果

Figure  2.  Isosurfaces of laser intensity distribution 3D phantoms and the corresponding reconstructions: (a−b) laser intensity with Gaussian distribution, kept unhanged along the propagation direction of the beam; (c−d) laser intensity with Gaussian distribution, changed along the propagation direction of the beam; (e−f) laser intensity with non-Gaussian distribution, unhanged along the propagation direction of the beam; (a), (c) and (e) are original phantoms and (b), (d) and (f) are the corresponding reconstructions

下载: 全尺寸图片 幻灯片

为了定量评估重建精度，采用公式（5）计算重建误差E ^[19]：

其中，E_g代表phantom，E_r代表重建结果。经计算，激光光强沿着光束传播方向无变化的高斯及非高斯分布的phantom重建误差分别为6.53%和5.87%,而沿着光束传播方向光强分布有变化的phantom重建结果为7.02%。在上述的仿真计算中，非高斯分布的重建精度略高于高斯分布，但并不能认为所有的非高斯分布phantom均能得到较高的重建精度，因为重建结果的准确性与模型的平滑度有关。值得注意的是，本工作旨在开发可重建真实实验情况中的激光光强三维分布算法，因此，仿真计算所采用的重建算法中并未加入激光分布的轴对称性假设。就上述计算结果来看，所采用的重建算法对不同类型的phantom均能得到较好的重建精度，具有一定的鲁棒性。因此，本工作提出的激光能量分布测量技术在原理上可行，即成像模型原理上正确，且所采用的重建算法收敛。

3.2   随机噪声对重建精度的影响

在实际测量中，不同信噪比的数据会对重建精度造成一定影响。为了验证模型在实际工况下的鲁棒性，进行数值仿真时采用单一变量法对激光光强呈高斯分布的phantom在多角度下的投影加入不同程度的噪声（2%、4%、6%、8%以及10%的随机噪声），并计算在不同信噪比下的重建误差。随机噪声的模型以及其数学公式为：

其中，$\vec{{{p}}}$即为投影，而$\vec{{{p}}}'$为添加随机噪声之后的投影。rnd为0−1之间的随机数，Noiselevel取值为2%~10%。

在重建区域离散程度相同的情况下，不同信噪比下激光光强三维分布phantom重建结果的等值面示意图及任一截面上激光光强分布示意图如图3（彩图见期刊电子版）所示。其中，（a−b）、（c−d）、（e−f）、（g−h）及（i−j）分别为施加2%、4%、6%、8%及10%噪声时的重建结果。同样地，为了定量评估重建精度，仍采用公式（4）计算重建误差。不同信噪比下的重建误差如图4中曲线所示。从图4可知，在施加10%以内的噪声时，重建误差均低于8.5%。

图  3  不同信噪比下激光光强分布三维重建等值面结果及中间截面激光光强二维分布示意图：（a−b）、（c−d）、（e−f）、（g−h）和（i−j）分别为施加2%, 4%, 6%, 8%和10%的随机噪声

Figure  3.  Isosurfaces of laser intensity distribution 3D reconstructions under different signal-to-noise ratios and the corresponding 2D slices for the middle section: (a−b), (c−d), (e−f), (g−h) and (i−j) are corresponding to adding 2%, 4%, 6%, 8% and 10% random noise, respectively

下载: 全尺寸图片 幻灯片

图  4  不同信噪比下的重建误差曲线

Figure  4.  Curve of the reconstruction error corresponding to different noise levels

下载: 全尺寸图片 幻灯片

4.   实验研究

4.1   实验方案

搭建一套基于层析成像的体激光诱导荧光信号探测系统，用于测量激光光强信号的三维分布。实验装置图如图5所示，通过将激光器固定在垂直于水平面的导轨上，使得激光束的传播方向为自上而下并垂直于水平面。比色皿中盛放有若丹明-乙醇溶液，当激光束（中心波长为532 nm）穿过若丹明-乙醇溶液时，会产生体激光诱导荧光信号。采用定制的一分多光纤束（购自南京春晖科技实业有限公司）同时采集7个角度的体激光诱导荧光信号。该光纤束有7个输入端，一个输出端。每个输入端近似布置于围绕比色皿所成的半圆周内，每个输入端前配有一个镜头（Nikon, f=50 mm, f/1.8），用于采集信号。光纤束输出端的信号经带通滤光片（中心波长为600 nm，FWHM=25 nm）传输至相机（Photron FASTCAM Mini AX100），相机前配有一个镜头（Nikon, f=55 mm, f/2.8）。经过滤光片滤除的信号可以有效滤去激光激发波长的影响。最后，将采集到的多角度体激光诱导荧光信号进行三维重建，可间接得到激光光强三维分布。

图  5  实验装置示意图

Figure  5.  Illustration of the experiment setup

下载: 全尺寸图片 幻灯片

各角度的位置关系如图6所示。图6（a）为各个角度距离重建区域的距离，图6（b）为在垂直于激光束的平面内绕激光束的各个角度之间的相对关系。7个角度分布在垂直于激光光束的平面近半圆周内，各个角度距离重建区域的距离约为160 mm，且7个角度的角度覆盖范围约为150°。

图  6  各角度位置示意图

Figure  6.  Schematic diagram of the position for different views

下载: 全尺寸图片 幻灯片

4.2   交替迭代三维重建

结合上述的层析成像技术原理以及建立的数学模型，需要考虑若丹明-乙醇溶液对体激光诱导荧光信号的吸收作用。结合比尔-郎伯定律^[24]，公式（2）可写为：

其中：α为混合均匀后若丹明-乙醇溶液的吸收率，可以通过测量激光束穿过溶液前后的功率变化，使用比尔-朗伯定律计算得到。

当同时记录不同角度待测区域内的体激光诱导荧光信号时，可获得多组线性方程，表示为向量形式：

由于方程非线性，本工作基于上述公式（4）代数重建法，采用一种交替迭代重建算法^[25]对体激光诱导荧光信号场进行求解：

其中：$\vec{{{f}}}({x}_{i},{y}_{i},{{\textit{z}}}_{i})$代表迭代过程中求得的体素$\left(x_{i}, y_{i}, {\textit{z}}_{i}\right)$的激光诱导荧光信号强度，上标k代表第k次迭代。λ_ART 为松弛因子，可控制迭代的收敛速度与收敛性, A为对应于像素$\left(x_{p}, y_{p}\right)$的权重系数。因此，在每一次迭代过程中都考虑了真实投影值与重建值间的差值。最终，通过设置迭代终止条件：迭代残差小于设定值或迭代次数大于设定值即可终止迭代。

4.3   结果与讨论

通过上述方法重建出的激光光强三维分布如图7（彩图见期刊电子版）所示。其中，图7（a）为三维结果示意图，图7（b）为激光束某中间截面的二维光强分布。

图  7  （a）三维激光光强重建结果示意图； （b）某中间截面激光光强二维分布示意图

Figure  7.  (a) 3D reconstruction results of the laser light intensity; (b) 2D center slice of the 3D reconstruction

下载: 全尺寸图片 幻灯片

没有验证方法可对测量误差进行直接评估，但是可通过计算实验测得的激光光强分布的二维投影数据与重建之后的反投影数据之间的相关性来证明重建算法的收敛性。与此同时，取实验测得的7个角度投影中的6个角度数据进行三维重建，计算余下一角度的重建反投影数据与实验测得的真实投影数据之间的相关性，由此间接验证重建方法的正确性^[26]。本工作中，考虑到激光直径较小，实验时在其周围摆放更多镜头较难实现。此外，已有学者对三维重建所需的角度数对重建结果的影响进行了较为详细的研究^[27]，研究结果表明，当重建角度大于等于5个时得到的重建结果在可接受范围之内。因此，本实验采用7个角度同时采集光强信号。其中，在进行重建时仅采用6个角度，利用余下一角度的重建反投影以及投影数据间的相关性来间接证明此重建方法的可行性。图8（a）（彩图见期刊电子版）为实验测得的7个角度的光强信号；图8（b）（彩图见期刊电子版）是未用于计算系数矩阵的某角度，其投影信号以及重建反投影信号的对比结果。采用Matlab软件工具箱中的函数 corr2计算投影以及重建反投影之间的相关系数，其计算原理如下：

图  8  投影与反投影对比示意图。（a）实验测得7个角度激光光强投影； （b）某角度投影与反投影（其中，投影和反投影相关系数为0.9802）

Figure  8.  Contrast of projection and re-projection. (a) Projections from seven views obtained experimentally；(b) projection and re-projection from a specific view (r=0.980 2)

下载: 全尺寸图片 幻灯片

其中，A和B分别代表投影和重建之后的反投影。经计算，该角度的投影以及反投影之间的相关性系数为0.9802，可以间接验证该方法的可行性。

5.   结　论

本文提出了一种新的基于层析成像技术的激光光强分布测量方法。该方法继承了层析成像技术的优点，同时借由内窥技术解决了传统层析成像技术设备造价高昂的问题。

通过数值仿真计算，对不同激光光强分布模型以及不同随机噪声等级时的重建精度进行评估。经计算，采用激光光强沿着光束传播方向无变化的高斯及非高斯分布的模型以及沿着光束传播方向光强分布有变化的模型，其重建误差分别为6.53%，5.87%和7.02%。在施加10%以内随机噪声时，重建误差均低于8.5%。即验证了重建模型的正确性以及重建算法的收敛性。此外，通过本工作提出的层析测量系统同时采集多个角度比色皿中体激光诱导荧光信号，并对其进行吸收矫正的三维重建，可间接获得激光光强分布。为了定量评估测量精度，在进行重建时仅采用其中6个角度，利用余下一角度的重建反投影以及投影数据间的相关性来间接证明此重建方法的可行性。经计算，该角度投影以及反投影之间的相关性系数为0.9802，可间接验证该方法的可行性。

在工业加工、激光传输及目标损伤效应试验等激光应用领域中均需要对激光能量分布参数进行准确测量。可以预见，本工作提出的激光光强三维测量方案在激光应用领域具有广泛的前景。