光声层析重建飞秒光丝二维横向图像仿真研究

准确获取飞秒激光成丝横截面图像及其沉积能量空间分布信息，对于成丝动力机制研究和促进诸多基于光丝的实际大气应用发展具有重要意义。本文基于热传导方程和波动方程构成的光声信号前向仿真模型，理论模拟了利用环阵式光声层析系统接收飞秒激光在空气介质中成丝诱导产生超声脉冲信号的过程；然后，利用延迟叠加算法对飞秒激光大气传输成丝沉积能量横向分布图像进行了反向重建，并分析了测量系统中关键器件超声换能器的中心频率、带宽、表面尺寸和探测表面灵敏度等性能参数对光丝沉积能量分布图像重建结果的影响。结果表明，单丝诱导产生的声压脉冲信号频谱为单峰结构，而多丝声压脉冲信号频谱为多峰结构；相比于单丝图像重建，多丝图像重建受“孔径效应”影响更显著；换能器的性能参数对光丝图像的重建效果有显著的影响，换能器的带宽越大、表面直径越小，以及表面灵敏度系数越大，越有利于光丝沉积能量分布图像重建效果的提升。该研究结果可为真实大气条件下飞秒激光传输成丝沉积能量空间分布的实验测量提供一定的理论支撑。

Abstract

Objective

Filament refers to a plasma channel with high laser intensity and high plasma density formed by the propagation of intense femtosecond laser pulses in a transparent medium. Several literatures have shown that the cross-section image of an optical filament at a specific z usually contains abundant structural information such as filament diameter, length, and energy distribution, which is of great significance for the visualization study of the dynamic process of filament formation. Moreover, accurate acquisition of the spatial structure and energy deposition distribution of femtosecond optical filaments are also of great significance for the development of filamentation-based atmospheric applications. Nevertheless, it is also the inherent parameter most difficult to measure directly. To solve the problem, we introduce a new medical imaging method named photoacoustic tomography (PAT) for optical filament cross-section imaging. The feasibility of reconstructing monofilament and multifilament images by photoacoustic tomography is verified theoretically. Moreover, we also study the influence of the performance parameters of the ultrasonic transducers on the optical filament image reconstruction.

Methods

We adopt a forward simulation model based on the photoacoustic wave equation to simulate the acquisition process of ultrasonic signals induced by optical filaments in air. A circular-scanning-based PAT system is considered to obtain the cross-section image of the laser filament. To simplify the problem, we assume that the initial heat source distribution of the optical filament satisfies the Gaussian distribution form, which can represent both the small high-energy core of the optical filament and its weak background energy region with a larger range. Based on experimental measurements, the initial maximum energy deposition density is assumed to be in the order of 10 mJ/cm³, and the diameter of the heat source is assumed to be in the order of 100 μm. The simulated time series of the acoustic signal is then applied to reconstruct the transverse distribution of femtosecond laser filaments with delay and sum (DAS) algorithm. Moreover, we also analyze the influence of performance parameters of ultrasonic transducers such as center frequency, bandwidth, surface size, and detection surface sensitivity on the reconstruction of filament cross-sectional images. The back-projection amplitude distribution profile along the y-axis is leveraged to compare the effect of image reconstruction.

Results and Discussions

According to the time series of ultrasound signals generated by monofilaments and multifilaments recorded at different detection distances, the frequency of monofilament and multifilament induced by femtosecond laser with multi-millijoule pulse energy is mainly concentrated within 4 MHz (Fig. 2). The signal spectrum of monofilament is single-peak structure, while the acoustic signal spectrum of multifilament is multi-peak structure (Fig. 2). The amplitude value of sound pressure signal decreases rapidly due to the attenuation of air. As the center of the optical filament deviates further from the scanning center, the cross-section image of the optical filament reconstructed by the back-projection (BP) algorithm and the DAS algorithm appears an obvious "elongated" phenomenon in the tangential direction (y-axis), which is the so-called "finite aperture effect" (Fig. 3). For monofilaments, the maximum energy amplitude decreases significantly with the increase in the center frequency of the transducer, which may be related to the filtering out of more low-frequency signals (Fig. 4). The same method is adopted to reconstruct the image of multifilament. It is found that the reconstructed multifilament image appears serious deformation with the multifilament center position deviating from the scanning center (Fig. 5). When x₀=1.0 mm, the two monofilaments near the scanning origin side can still be distinguished, whereas the two monofilaments near the transducer side are fused and cannot be distinguished. Therefore, the secondary filaments around the multiple filaments are more susceptible to the "aperture effect" and the fuzzy deformation occurs. The fuzzy deformation effect will be more obvious when the distance becomes larger from the scanning center or the distance becomes smaller from the surface of the transducer. Therefore, compared with monofilament reconstruction, multi-filament image reconstruction is more affected by the "aperture effect". Especially, the blur deformation of the surrounding sub-filaments is more likely. In summary, the characteristics of the transducer have an obvious influence on the reconstruction of monofilament and multifilament cross-sectional images. A larger bandwidth of the transducer will cause a smaller surface diameter, a larger surface sensitivity parameter, and a better reconstruction quality of monofilament and multifilament images. The influence of the center frequency of the transducer on the optical fiber image reconstruction is very complicated. Therefore, it is necessary to select the transducer with the appropriate center frequency combined with the spectrum analysis of the acoustic signal in the actual measurement.

Conclusions

We utilize a novel medical imaging method named PAT to reconstruct cross-section images of femtosecond laser filament formed in an air medium. The results show that the acoustic signal induced by a single filament has a single-peak structure, while that induced by a multifilament has a multi-peak structure. The performance parameters of the transducer have an obvious influence on the reconstruction results. A larger bandwidth of the transducer will lead to a smaller surface diameter, a larger surface sensitivity coefficient, and a better reconstruction effect of energy deposition distribution of optical filament. Compared with monofilament, the reconstruction of the multifilament image is more susceptible to the "finite aperture effect". Our study can provide some theoretical support for the experimental measurement of the spatial deposited energy distribution of femtosecond laser filament transmission under real atmospheric conditions.

1　引言

高功率的飞秒激光脉冲在大气中传输时，由于克尔自聚焦、衍射和等离子体散焦效应之间的动态平衡，可以形成长度很长的等离子体通道，即光丝^［1］。飞秒激光成丝包含丰富的物理过程，并在远距离大气探测、引导放电、产生有机物和诱导水汽凝结等诸多大气研究领域中展现出潜在的应用价值，逐渐成为非线性光学与大气科学交叉研究的前沿热点问题之一^［2-6］。

空气中形成的等离子体光丝，其重要特性之一就是它能够几乎在瞬间沉积相当一部分脉冲能量到传输介质^［7］。由于光场电离和空气分子的非共振转动拉曼激发等作用是造成光丝能量沉积的主要物理过程^［7］，沉积能量的分布与光丝内的光场强度和自由电子密度紧密相关。因此，通过获取光丝沉积能量空间分布，可以间接了解自由电子的分布和光丝的空间结构，这对于成丝动力学过程的“可视化”研究具有重要意义^［8］。同时，这部分沉积的脉冲能量，最终以热量的形式被介质吸收，能够激发形成压力脉冲波以及更长时间尺度的热传导过程，相关物理过程也被认为是诱导产生气流扰动、冰晶繁生和清理光学通道等现象的关键^［9-11］。然而，由于光丝内部强度极高，使得目前对于光丝能量沉积密度空间分布的直接测量还非常困难。现有诊断光丝的方法大多是基于成丝过程中伴随产生的光、声和热等信号间接实现的^［12］。其中，最为常用的是利用微麦克风测量声信号强度来间接估算能量沉积效率^［7，13］。该方法实施简单，测量过程中只需要沿着光丝传输方向移动高频微麦克风，获取不同位置的声信号强度，然后根据热能与声压之间的关系 $ε = (5 / 2) p$ 就能间接估算光丝线性能量沉积率。然而，这种基于微麦克风的测量方法，空间精度太低，无法获得与光丝内部结构有关的参数^［14］。

光声层析（PAT）是一种新兴的非侵入式医学成像技术，它利用光声效应对介质（如生物组织）中沉积的光能量分布进行成像^［15］。利用PAT技术重建单丝横截面图像的可行性，在理论和实验上都已经得到了充分验证^{［14，16-17］}，并且还逐步向着三维重构的趋势发展^［18］，是一种非常有发展前景的定量获取光丝结构信息的探测方式。在PAT实验测量系统中，超声换能器是接收光声信号的关键元件，对最终重建图像的质量有重要的影响。实际实验中所用的换能器通常具有有限的尺寸，一般在3~10 mm之间，并且换能器的带宽也是有限的^［19］。Xu等^［20］通过推导，发现换能器的有限孔径和带宽对PAT成像的切向分辨率有明显的影响。另外，在实际测量中，所用换能器的中心频率和表面灵敏度也是有限的，接收的声信号还会存在各种噪声。这些因素综合起来对光丝图像重建的影响如何，非常值得进一步深入研究^［21］。

为此，本文通过综合考虑换能器的中心频率、带宽、尺寸和灵敏度等性能参数，理论模拟了利用特定环阵式PAT系统探测单丝和多丝声信号的过程，并验证了利用延迟叠加算法重建出光丝能量沉积分布横截面图像的可行性。该研究结果可为进一步重建光丝沉积能量三维空间分布奠定理论基础，对诸多基于光丝的大气应用研究具有较好的参考价值。

2　研究方法

PAT成像大致包括两个过程：前向过程和后向过程。前向过程指的是高功率短脉冲激光照射吸收体过程中，由于热弹性膨胀激发出宽带超声波信号，并被超声换能器接收的过程。后向过程是指利用一定图像重建算法通过求解逆问题，根据收集到的光声信号，重建目标物光学吸收能量分布或者初始声压分布的过程^［15］。

2.1　光丝诱导超声信号前向仿真模型

光丝诱导超声信号前向仿真模型的建立，需要结合光丝诱导产生声信号所经历物理过程的时间尺度考虑。在飞秒激光脉冲通过传输介质后，先激发电子密度的相干振荡过程（~100 ps），形成等离子体尾迹；产生的等离子体及其吸收的部分脉冲能量在等离子体复合之后，最终转换成气体的热能（~10 ns）；随后，为达到机械平衡，会产生压力脉冲波，激发形成声信号（~1 μs）；在压力平衡后，发生热传导过程，可持续至1 ms以上量级^{［7，13，22］}。由此可见，飞秒光丝激发形成声信号的过程，本质上是一个热致声过程。近似地，光丝诱导产生的初始压力与光丝沉积的热量成正比^［13］。因此，通过重建光丝诱导初始声压强度分布，可以间接得到相应的沉积能量分布^［23］。根据对相关物理过程的时间尺度的分析结果，在考虑热扩散，忽略介质黏性的情况下，前向过程可以用下列方程组^［24］进行描述：

\{\begin{matrix} \frac{\partial}{\partial t} (T_{1} - \frac{γ - 1}{γ α} p_{1}) = χ Δ T_{1} + \frac{H}{ρ_{0} c_{p}} \\ (Δ - \frac{1}{v_{s}^{2}} \frac{\partial^{2}}{\partial t^{2}}) p_{1} = - ρ_{0} β \frac{\partial^{2} T_{1}}{\partial t^{2}} \end{matrix}

，（1）

式中：p₁和T₁分别为气压和温度的扰动量；H为初始热源分布；ρ₀为介质密度；v_s为绝热声速； $γ = c_{p} / c_{v}$ 为比热容；c_p为定压比热容；c_v为空气定容比热容； $χ$ 为热扩散率；β为热膨胀系数；Δ表示拉普拉斯梯度算子。

在忽略热扩散、介质的黏度，以及沿z方向激光辐射强度衰减的情况下，联立式（1）^［14］可得

\frac{1}{r} \frac{\partial}{\partial r} (r \frac{\partial p}{\partial r}) - \frac{1}{v_{s}^{2}} \frac{\partial^{2} p}{\partial t^{2}} = - \frac{β}{c_{p}} \frac{\partial \tilde{H}}{\partial t}

，（2）

式中： $p (\vec{r}, t)$ 表示位于原点O的 $\vec{r}$ 位置上的换能器接收到的声压信号时间序列； $\tilde{H} (\vec{r}, z, t) = H (\vec{r}) δ (t)$ 为激光脉冲作用函数， $δ (t)$ 为Dirichlet函数。为了得到不同位置和不同时刻的声压脉冲信号，需要以初始光丝能量吸收分布为输入条件。根据以往的研究，单丝横截面分布通常由单个能量核心以及周围相对较大区域的弱背景能量区构成，而多丝横截面表现为多个随机分布的“热点”^［25］。因此，为了简化问题，可假设单条光丝沉积的脉冲能量被介质以热量形式吸收后产生的初始热源满足高斯分布^［22］：

H (r) = H_{0} e x p [- {(r / w_{0})}^{2}]

，（3）

式中：H₀表示单位体积内激光脉冲沉积的最大能量；w₀为光丝热源分布高斯半径；r为径向坐标位置。空气介质中初始气压分布与沉积脉冲能量之间近似满足如下关系^［16］：

p (r, t = 0) = (\frac{β v_{s}^{2}}{c_{v}}) H (r)

。（4）

将式（3）和式（4）代入式（2）并进行傅里叶变换，可得

\tilde{p} (\vec{r}, f) = \frac{π^{2} w^{2} H_{0} f β}{2 c_{p}} H_{0}^{(2)} (\frac{2 π f}{v_{s}} r) e x p [- {(\frac{π f w}{v_{s}})}^{2}]

，（5）

式中： $\tilde{p} (\vec{r}, f)$ 为 $p (\vec{r}, t)$ 的傅里叶变换量；f为角频率； $H_{0}^{(2)} (u)$ 为Hankel函数。当接收信号位置 $r ≫ w_{0}$ 时，进一步可以得到轴对称坐标系下的声压信号 $p (r, t)$ ^［26］为

p (r, t) = \frac{H_{0} w_{0}^{\frac{1}{2}}}{2 \sqrt[]{2 π r}} \frac{β v^{2}}{c_{p}} F (\frac{r - v_{s} t}{w_{0}})

，（6）

F (ζ) = [Γ (\frac{3}{4}) {}_{1}F F_{1} (- \frac{1}{4}; \frac{1}{2}; ζ^{2}) + 2 ζ Γ (\frac{5}{4}) {}_{1}F F_{1} (\frac{1}{4}; \frac{3}{2}; ζ^{2})] e x p (- ζ^{2})

，（7）

式中： $Γ (\cdot)$ 为Gamma函数； ${}_{1} F_{1} (a; b; ζ)$ 为合流超几何函数。

在后续的仿真中，主要的参数取值为v_s=330 m/s，c_p=1.002 kJ/（kg·K），c_v=0.717 kJ/（kg·K），β=1/273.15 ℃^-1。H₀和w₀的取值在后文中说明。另外，由于空气介质的黏性，超声信号在空气中传输时会发生明显的衰减。为了考虑声压信号的衰减作用，文献［24］中提出了一种等效方法，即将热源半径w₀等效为 $\sqrt[]{w_{0}^{2} + (v_{s} {/ π)}^{2} α r}$ ，其中r为传播距离， $α \approx 1.85 \times 10^{- 11}$ m^-1Hz^-2，本文也采用该方法进行计算。

2.2　图像重建方法

光声图像重建是一个典型的逆问题，即由光声信号p计算出介质电磁吸收分布A的问题。在诸多的典型PAT实验系统中，2D环形阵列扫描可以覆盖较为完整的目标视角，常用于吸收体断层面的成像，其示意图如图1（a）所示。脉冲激光照射在目标物上激发产生超声信号，换能器围绕目标物沿着环形轨迹扫描或者使用多个换能器组成环形阵列，在不同角度位置上接收超声信号。然后，根据不同位置换能器接收到的信号，利用一定的图像算法重建出初始的光吸收分布。对于环阵扫描，最常用的图像重建算法是延迟叠加（DAS）算法。该算法的实质是已知所有换能器采集的信号值，根据图像重建区域像素点与换能器之间的距离大小，求出某一像素与每个换能器之间的时间延迟，然后由相应时间延迟的换能器信号之和得出该点像素值。同时，通过融合立体角权重因子，可有效补偿检测视点的变化。加权形式的延迟叠加算法计算公式^［27］可以表示为

p_{0}^{(b)} (\vec{r}) = \int_{Ω_{0}} b ({\vec{R}}_{i}, t = \frac{|\vec{r} - {\vec{R}}_{i}|}{v_{s}}) \frac{d Ω_{0}}{Ω_{0}}

，（8）

式中： $\vec{r}$ 为扫描中心点到重建像素点的距离； ${\vec{R}}_{i}$ 为扫描中心点到第i个换能器的距离； $p_{0}^{(b)} (\vec{r})$ 为 $\vec{r}$ 位置处重建得到的声压值； $b ({\vec{R}}_{i}, t) = 2 p ({\vec{R}}_{i}, t) - 2 v_{s} t \partial p ({\vec{R}}_{i}, t) / \partial t$ 表示反投影项； $p ({\vec{R}}_{i}, t)$ 表示第i个换能器接收到的声信号，根据式（6）计算得到； $Ω_{0}$ 为整个测量面S₀相对于其内部 $\vec{r}$ 位置重建点的立体角。对于平面， $Ω_{0} = 2 π$ ；对于球形和圆柱体， $Ω_{0} = 4 π$ ； $d Ω_{0} / Ω_{0}$ 为立体角加权因子， $d Ω_{0}$ 为换能器表面检测单元dS₀相对于 $\vec{r}$ 位置重建点的立体角：

d Ω_{0} = \frac{d S_{0}}{{|\vec{r} - {\vec{R}}_{i}|}^{2}} [{\vec{n}}_{0} \cdot \frac{(\vec{r} - {\vec{R}}_{i})}{|\vec{r} - {\vec{R}}_{i}|}]

。（9）

在DAS算法中，为了综合考虑换能器的中心频率、带宽、尺寸和灵敏度等性能参数，如图1（b）所示，对于一个表面长度为L的换能器，将其表面均匀地划分成N个探测单元，然后将这N个探测单元接收的信号合成一个信号后进行反投影^［28］。通常情况下，换能器表面的接收灵敏度并不均匀，一般由中心向边缘逐渐减小。因此，本文采用高斯函数来表示换能器表面接收声压的灵敏度情况，相应的灵敏度系数^［28］可以表示为

A ({\vec{R}}_{i}^{k}) = e x p (- \frac{{|{\vec{R}}_{i}^{k} - {\vec{R}}_{i}^{}|}^{2}}{2 σ^{2}})

，（10）

式中： $σ$ 为换能器表面灵敏度系数； ${\vec{R}}_{i}^{k}$ 扫描中心点到第i个换能器表面第k个探测单元的向量，探测单元距离换能器表面中心点越远，相应灵敏度系数就越小。

图 1. 环阵PAT仿真模型设置。（a）环形阵列探测示意图；（b）换能器表面探测单元划分示意图

Fig. 1. Simulation model setup of ring array PAT system. (a) Schematic diagram of ring array detection; (b) schematic diagram of division of transducer surface detection unit

下载图片查看所有图片

作为对比，本文还采用传统的反投影（BP）算法对光丝图像进行重建。所用到的BP算法有两种：一种是理想点换能器模型（记为BPoint）；另一种是无限大平面换能器模型（记为BPinf）。BPinf算法的反投影线是一组与探测器平面平行的直线，相应的反投影值^［19］为

I (x, y) = \sum_{i = 1}^{n} S_{i} [(R - x c o s φ_{i} - y s i n φ_{i}) / v_{s}]

，（11）

式中： $S_{i} (t)$ 表示第i个换能器中心接收到的超声信号； $φ_{i}$ 表示该换能器与x轴的夹角，第一个换能器设定位于x轴正向上； $I (x, y)$ 表示重建区域 $(x, y)$ 像素点上的反投影幅度值。

3　分析与讨论

3.1　光丝声信号特征

实验测量结果表明，由Ti：Sapphire飞秒激光器（805 nm，55 fs，10 Hz）产生的入射能量E_in约为2~3 mJ的单脉冲在空气介质中传输成丝时，形成的单丝内部能量沉积密度峰值H₀=（50±10）mJ/cm³，光丝能量沉积分布半径w₀=（60±10）μm^［29］。如果在光束传输路径上增加掩膜板，相应的能量沉积密度变为H₀=（15±5）mJ/cm³，光丝能量沉积分布半径w₀=（250±50）μm。因此，本文假设由初始单丝沉积的能量分布满足式（3）所示的高斯分布，取H₀=60 mJ/cm³，w₀=100 μm/（ $2 \sqrt[]{l n 2}$ ），即半峰全宽（FWHM）为100 μm，结果如图2（a）所示。根据式（4）计算可得，该单丝结构诱导产生的初始气压扰动最大幅度值在20 kPa以上，该量级符合激光等离子体激波压强取值^［30］。从图2（a）还可以看出，高斯分布在一定程度上既可以表示光丝核心高强度区域，也能表示光丝周围的弱背景能量池^［31］。假设模拟的多丝结构由5条单丝按照正方形排列组合而成。其中，能量较大的主丝［H₀=24 mJ/cm³，w₀=60 μm/（ $2 \sqrt[]{l n 2}$ ）］位于正方形中心，4条能量较小的副丝［H₀=12 mJ/cm³，w₀=60 μm/（ $2 \sqrt[]{l n 2}$ ）］位于正方形4个角上，副丝与主丝之间的距离为 $150 \sqrt[]{2}$ μm，如图2（d）所示。多丝结构诱导产生的初始气压扰动最大幅度值在9 kPa以上。

图 2. 飞秒光丝光声仿真结果。（a）初始单丝沉积能量分布横截面；（b）不同接收位置处的声压信号；（c）声压信号的频谱分布；（d）初始多丝沉积能量分布横截面；（e）不同接收位置处的声压信号；（f）声压信号的频谱分布

Fig. 2. Photoacoustic simulation results of femtosecond laser filament. (a) Initial cross-sectional energy deposition distribution of single filament; (b) acoustic pressure signals recorded at different positions; (c) frequency spectrum of acoustic pressure signals; (d) initial cross-sectional energy deposition distribution of multiple filaments; (e) acoustic pressure signals recorded at different positions; (f) frequency spectrum of acoustic pressure signals

下载图片查看所有图片

图2（b）和2（e）分别为不同探测距离上记录的单丝和多丝产生的声信号时间序列。图2（c）和2（f）分别为声信号相对应的频谱分布。从图2（b）中可见，单丝产生的声信号具有典型的“N”字型结构，对应于一个陡峭的正压脉冲以及紧随其后的一个负压脉冲，这与文献［24］中实际测量结果一致。随着传输距离的增大，声压信号最大幅度值相比于初始20 kPa值已经明显减小，但波形几乎没有变化。利用换能器在5 mm位置处接收声信号，最大声压幅度值就已经减小到了1 kPa以下。多丝产生的声信号具有明显的多个峰值结构［图2（e）］。这种多峰结构的形成，与多丝内部结构中不同强度中心与换能器距离不同相关。从图2（c）和2（f）所示声信号频谱来看，单丝和多丝产生的声信号频率主要集中在4 MHz以内，单丝声信号频谱为单峰结构，而多丝声信号频谱为多峰结构。距离越远，声信号高频部分衰减得越明显，因此空气的衰减作用相当于产生了低通滤波效应，滤除了时域信号中的高频成分，使得信号频谱峰值向低频偏移。声压信号幅度值的迅速减小就与高频信号的衰减相关。但是由Kramers-Kronig关系可知，空气中的声速可以认为保持不变^［32］。单丝诱导产生的声脉冲信号频谱峰值约为0.8 MHz，比文献［24］报道的实际测量结果1.16 MHz要略小。这主要与本文探测距离（5 mm）比文献中探测距离（~3 mm）更大，声信号传输过程中衰减较大有关。对比来看，在考虑声信号衰减的情况下，多丝声信号最大峰值频率（约2.0 MHz）比单丝峰值频率（约0.8 MHz）更高，单丝声信号高频部分的衰减更为明显。由于高频段的光声信号能突出物体的细微结构，因此为了能够更有效地接收光丝诱导形成的声信号和重建出光丝图像，需要使用工作频率达到兆赫兹量级的空气耦合式换能器进行测量。

3.2　单丝图像重建

在仿真得到单丝和多丝产生的声信号以后，就可以利用BP和DAS两种图像重建算法对光丝横向沉积能量分布图像进行重建。具体实施过程中，先重建得到初始气压分布，然后利用式（4）转换成吸收能量分布。为了更接近实际的测量过程，进一步假设数据采集卡的采样频率为50 MHz，环阵探测所用的换能器数量为64个，也可以等价为单个换能器旋转扫描，并在64个位置上记录采集到的声信号。单个换能器表面长度为3 mm，带宽和中心频率分别为70%和2.0 MHz，并且对声信号添加40 dB的随机噪声。其中，BPoint和BPinf重建使用的是换能器中心点对应的声信号时间序列。而在实施DAS算法时，是先将换能器表面划分成51个探测点，然后根据光丝热源中心到每个探测点的距离计算得到每个探测点接收到的声信号，最后再将这51个探测点对应的声信号进行叠加，得到合成信号后再重建。由于在BP算法中，反投影的是换能器中心一个点接收的声信号，因此为了和DAS算法进行对比，将理想点换能器模型重建结果放大了51倍。在合成DAS信号过程中，还假设换能器表面接收灵敏度满足高斯分布，也就是换能器中心接收信号能力强，边缘接收能力弱，这也更符合实际测量情况。高斯分布的最大FWHM为σ=0.5 mm。环阵扫描半径为5 mm，环阵内部重建区域x方向总长度为2 mm，y方向总长度为2 mm，网格点为1024×1024，分辨率约为2.0 μm。

图3所示为当光丝中心分别位于x轴上不同位置时，利用BP和DAS算法重建得到的单丝沉积能量分布横截面图像。其中x₀表示光丝中心距离扫描原点O的距离。第一列、第二列和第三列分别为BPoint、BPinf和DAS图像算法重建得到的光丝横截面图像。第四列所示为沿重建截面图像y轴方向的能量密度幅值分布。从重建的横截面分布图来看，当光丝核心位于扫描原点O时［图3（a）］，BP和DAS两种图像算法均能很好地重建出的横截面分布轮廓。尤其是能重建出光丝核心周围的弱“能量池”，这是其他干涉法、电导法等大多数间接测量方法无法得到的^［16］。光丝重建图像周围还存在能量密度为负值（也就是气压为负值）的区域，这与滤除了低频信号有关^［16］。相比于初始假设的H₀=60 mJ/cm³，3种算法重建得到的初始能量密度最大值均更小。造成这种差异的原因一方面在于虽然本文在前向仿真中考虑了超声信号衰减，但是在后向反投影过程中没有考虑衰减补偿；另一方面在于本文使用的是具有一定带宽和中心频率的超声换能器，接收信号过程中起到了一定的滤波作用，部分高频信号被滤除。最终，使得重建得到的能量密度值大幅度减小。对比来看，BPinf算法重建得到的能量密度最大值在三者中最大，更接近初始值，由此可知，使用更大探测面的换能器有利于提高重建幅度值。这可能与5 mm扫描距离位于换能器近场范围内，声场传播符合直线传播规律有关。在实际应用中，也通常采用探测面积较大的非聚焦平面超声换能器来解决信噪比低的问题^［33］。但是，探测面太大的换能器更容易受到“孔径效应”的影响^［27］。这可以从图3（b）和3（c）中看出。随着光丝偏离扫描中心位置越来越远，BPinf算法和DAS算法重建出的光丝横截面图像在切向方向（y方向）上，出现了明显的“被拉长”现象，这就是“孔径效应”。当光丝中心远离扫描中心时，重建效果变差。尤其是当x₀=2.0 mm时，BPinf算法重建图像基本上不能识别出单丝。从第四列重建廓线来看，还出现了多个峰值中心。这说明“孔径效应”对于探测面大的换能器表现得更加明显。对于不同位置上的光丝，BPoint算法基本上能够较好地重建出光丝轮廓信息，因此使用具有小尺寸探测面的换能器有利于重建出光丝的细微结构。总的来看，在进行光丝测量时，应该尽可能地将光丝置于扫描中心位置。

光丝中心分别位于x轴上不同位置时单丝沉积能量分布横截面图像重建结果。（a）x0=0 mm；（b）x0=1.0 mm；（c）x0=2.0 mm

图 3. 光丝中心分别位于x轴上不同位置时单丝沉积能量分布横截面图像重建结果。（a）x₀=0 mm；（b）x₀=1.0 mm；（c）x₀=2.0 mm

Fig. 3. Reconstructed cross-sectional energy deposition distribution of single filament when centers of filament are located at different positions on x-axis. (a) x₀=0 mm; (b) x₀=1.0 mm; (c) x₀=2.0 mm

下载图片查看所有图片

进一步利用DAS算法，综合分析了换能器的扫描半径D_r、换能器数量N、换能器的中心频率f_c、换能器带宽W、换能器探头直径L和表面灵敏度系数σ对单丝重建结果的影响，结果如图4所示。具体计算参数的设置，可以参见表1。从图4（a）可见，随着环阵扫描半径的增大，重建初始沉积能量的最大幅度值略有减小。这可能是由于随着扫描半径的增大，换能器的探测灵敏度下降所导致的^［34］。环阵换能器数量对单丝重建几乎没有影响［图4（b）］。换能器的中心频率对重建图像的幅度值有明显影响［图4（c）］。对于单丝，随着换能器中心频率增加，最大能量幅度值明显减小，这与较多低频信号被滤除有关。换能器的带宽、直径和表面灵敏度等特性参数对重建图像也有明显的影响［图4（d）、4（e）和4（f）］。换能器的带宽W越大、直径L越小和表面灵敏度参数σ越大，重建能量分布最大幅度值也越大。

图 4. 不同换能器相关因素对单丝沉积能量分布重建结果的影响。（a）扫描半径D_r；（b）换能器数量N；（c）换能器中心频率f_c；（d）换能器带宽W；（e）换能器直径L；（f）换能器表面灵敏度σ

Fig. 4. Influence of different transducer related factors on reconstruction of energy deposition distribution of single filament. (a) Scanning radius D_r; (b) number of transducers N; (c) center frequency f_c; (d) transducer bandwidth W; (e) transducer diameter L; (f) transducer surface sensitivity σ

下载图片查看所有图片

表 1. 不同参数下重建图像切向分辨率变化值

Table 1. Changes in tangential resolution of reconstructed image under different parameters

Factor	Symbol	Value	FWHM /μm	Other parameter
Scanning radius		5 mm	98.2	N=64，L=3 mm，f_c=2 MHz，W=70%，σ=1.0 mm
	D_r	8 mm	102.9
		10 mm	103.6
Number of transducers		32	97.4	D_r=5 mm，L=3 mm，f_c=2 MHz，W=70%，σ=1.0 mm
	N	64	98.2
		128	98.6
Transducer diameter		2 mm	100.6	D_r=5 mm，N=64，f_c=2 MHz，W=70%，σ=1.0 mm
	L	4 mm	98.4
		6 mm	98.7
Center frequency		1 MHz	170.2	D_r=5 mm，N=64，L=3 mm，W=70%，σ=1.0 mm
	f_c	2 MHz	98.2
		3 MHz	83.4
Transducer bandwidth		60%	107.0	D_r=5 mm，N=64，L=3 mm，f_c=2 MHz，σ=1.0 mm
	W	70%	98.2
		80%	96.4
Transducer surface sensitivity		0.5 mm	96.5	D_r=5 mm，N=64，L=3 mm，f_c=2 MHz，W=70%
	σ	1.5 mm	98.5
		3.0 mm	98.4

查看所有表

表1给出了重建光丝热源切向剖面的FWHM值，也就是切向分辨率值。总体来看，DAS算法重建出的单丝热源FWHM值与实际初始值（100 μm）较为接近。随着换能器数量的增加，切向分辨率值略有提升，更接近100 μm。换能器表面宽度越小，切向分辨率值越接近100 μm。换能器中心频率和带宽对切向分辨率值影响较大，换能器中心频率和带宽较小时，单丝重建结果切向分辨率值明显偏大，即图像出现明显的模糊^［27］。综合来看，为了提高重建图像质量，应该使用中心频率适中、带宽较大、直径较小和表面灵敏度参数较大的换能器。

3.3　多丝图像重建

利用BP和DAS两种算法对多丝图像进行重建，结果如图5所示。图中前三列分别是利用算法重建出的多丝沉积能量分布横截面图像，第四列为重建图像的对角线幅度值分布曲线。其中，采集卡的采样频率为50 MHz，环阵换能器总数为64个。单个换能器表面长度为3 mm，带宽和中心频率分别为70%和2.0 MHz，并且对声信号添加40 dB的随机噪声。同样地，在执行DAS算法重建时，将每个换能器表面划分成51个探测点，换能器表面接收灵敏度系数σ=0.5 mm。环阵扫描半径为5 mm，环阵内部重建区域x方向总长度为2 mm，y方向总长度为2 mm，网格点为1024×1024，分辨率约为2.0 μm。

光丝中心分别位于x轴上不同位置时多丝沉积能量分布横截面图像重建结果。（a）x0=0 mm；（b）x0=1.0 mm；（c）x0=2.0 mm

图 5. 光丝中心分别位于x轴上不同位置时多丝沉积能量分布横截面图像重建结果。（a）x₀=0 mm；（b）x₀=1.0 mm；（c）x₀=2.0 mm

Fig. 5. Reconstructed cross-sectional energy deposition distribution of multiple filaments when centers of filaments are located at different positions on x-axis. (a) x₀=0 mm; (b) x₀=1.0 mm; (c) x₀=2.0 mm

下载图片查看所有图片

从图5（a）可见，当多丝的主丝中心位于扫描中心点时，BP算法和DAS算法基本上都能够较为清晰地重建出主丝和周围副丝的能量分布轮廓。由于滤除了部分低频信号，主丝和副丝周围都出现了明显的负值区域。从图5（b）和5（c）可见，随着多丝中心位置偏离扫描中心，重建得到的多丝图像出现了严重的形变。当x₀=1.0 mm，靠近扫描原点一侧的两条单丝还能够分辨出来，但是靠近换能器一侧的两条单丝基本上融合在一起，无法区分。因此，多丝周围的副丝更容易受“孔径效应”的影响而发生模糊变形，并且越偏离扫描中心（或靠近换能器表面），这种模糊变形作用越明显。当x₀=2.0 mm，BPinf算法已经基本上区分不了主丝和副丝的位置，DAS算法重建得到的副丝也混淆在一起，所得幅度值没有太大差别。因此，在进行多丝图像光声探测时，一方面应该尽可能地将多丝中心置于扫描中心，另一方面还需要进一步从重建算法入手进行改进^［33］。

同样地，进一步利用DAS算法，研究了扫描半径D_r、换能器数量N、换能器的中心频率f_c、换能器带宽W、换能器探头直径L和表面灵敏度系数σ等参数对多丝图像重建结果的影响，结果如图6所示。仿真过程中，除了可变参数外，其他仿真参数均从D_r=5 mm、N=64、f_c=2.0 MHz、W=70%、L=3 mm和σ=1 mm这些参数中选取。比如在图6（a）中，D_r为可变参数，依次设置为5、10、15 mm，而其他参数则设置为N=64、f_c=2.0 MHz、W=70%、L=3 mm和σ=1 mm。从图6（a）可见，随着扫描半径的增大，沉积能量幅度值分布形状没有太大变化，但是主丝和副丝的最大幅度值逐渐增大，表明多丝探测需要较大的扫描半径。从图6（b）可见，换能器数量对多丝重建结果影响不大，这与单丝一致。从图6（c）可见，换能器中心频率也对多丝重建有明显的影响。当f_c=2.0 MHz时，多丝重建效果最好，其次是3 MHz，因此实际测量时需要结合声信号频谱分析选择具有合适中心频率的换能器，这样损失的频率成分较少，成像效果比较好。从图6（d）、6（e）和6（f）可见，换能器的带宽W越大、直径L越小和表面灵敏度参数σ越大，重建得到的沉积能量最大幅度值越大，也就越接近初始设置吸收能量分布，这与单丝重建的结论相一致。

图 6. 不同换能器相关因素对多丝沉积能量分布重建结果的影响。（a）扫描半径D_r；（b）换能器数量N；（c）换能器中心频率f_c；（d）换能器带宽W；（e）换能器直径L；（f）换能器表面灵敏度σ

Fig. 6. Influence of different transducer related factors on reconstruction of energy deposition distribution of multiple filaments. (a) Scanning radius D_r; (b) number of transducers N; (c) center frequency f_c; (d) transducer bandwidth W; (e) transducer diameter L; (f) transducer surface sensitivity σ

下载图片查看所有图片

4　结论

基于热传导方程和波动方程构成的光声信号前向仿真模型，分别研究了飞秒激光在大气中传输形成单丝和多丝结构诱导产生声信号的特征，并利用BP算法和DAS算法对单丝和多丝沉积能量横截面分布图像进行了重建。初步验证了利用环阵PAT重建单丝和多丝沉积能量横向分布图像的可行性，并得到以下结论：

1）在考虑声信号衰减的情况下，单丝诱导产生的声压信号为单峰结构，而多丝声信号为多峰结构；随着传播距离的增大，声信号高频部分衰减明显，单丝峰值频率（约0.8 MHz）要比多丝的最大峰值频率（约2.0 MHz）小；

2）PAT法具备重建多丝图像的能力，但是相比于单丝重建，多丝图像重建受“孔径效应”的影响更大，尤其是周围副丝更容易出现模糊变形；

3）换能器特性对单丝和多丝横截面图像的重建有明显影响，换能器的带宽越大、表面直径越小和表面灵敏度参数越大，越有利于提高单丝和多丝图像重建质量。换能器中心频率对光丝图像重建的影响较为复杂，实际测量时需要结合声信号频谱分析选择具有合适中心频率的换能器。结合本文研究结果来看，对数毫焦飞秒脉冲产生的单丝和多丝，当最大沉积脉冲能量密度在数十毫焦每立方厘米量级时，选择中心频率为2 MHz的超声换能器测量效果较好。

值得注意的是，DAS算法本质上也是一种BP算法，而BP方法是一种解析计算方法，虽然计算速度快，但缺乏有效获取成像区域定量信息的能力。在本文研究中，虽然DAS算法能够较好地重建出单丝和多丝的沉积能量横向分布轮廓信息，并且也具有一定的抗噪能力，但是最后重建出的光丝沉积能量最大幅度值与初始值存在较大差异。这种结果的影响因素是多方面的，主要与超声信号在空气中剧烈衰减，以及换能器的中心频率和带宽有限相关。后续可以使用基于模型的重建算法来综合考虑传输介质特性、声速衰减补偿和换能器特性等影响因素，从而更好地估计出光丝沉积能量分布。另外，开展相关的测量实验也将是下一步的工作重点。

参考文献

[1] Couairon A, Mysyrowicz A. Femtosecond filamentation in transparent media[J]. Physics Reports, 2007, 441(2/3/4): 47-189.

[2] 王铁军, 陈娜, 郭豪, 等. 飞秒强激光大气遥感新技术的原理和研究进展[J]. 激光与光电子学进展, 2022, 59(7): 070001.

Wang T J, Chen N, Guo H, et al. Principle and research progress of atmospheric remote sensing by intense femtosecond lasers[J]. Laser & Optoelectronics Progress, 2022, 59(7): 070001.

[3] 曾庆伟, 高太长, 刘磊, 等. 飞秒激光成丝诱导形成水凝物的机理研究进展[J]. 红外与激光工程, 2019, 48(4): 0406002.

Zeng Q W, Gao T C, Liu L, et al. Advances in mechanism research of femtosecond laser filamentation induced hydrometeors formation[J]. Infrared and Laser Engineering, 2019, 48(4): 0406002.

[4] Wolf J P. Short-pulse lasers for weather control[J]. Reports on Progress in Physics, 2018, 81(2): 026001.

[5] 冯志芳, 刘勋, 郝婷, 等. 面向遥感应用的飞秒激光脉冲超长距离传输的研究进展[J]. 中国激光, 2023, 50(7): 0708003.

Feng Z F, Liu X, Hao T, et al. Review on ultra-long distance propagation of femtosecond laser pulses for remote sensing applications[J]. Chinese Journal of Lasers, 2023, 50(7): 0708003.

[6] 程俊皓, 胡理想, 王铁军, 等. 飞秒激光诱导多丝产生与调控研究进展[J]. 中国激光, 2023, 50(14): 1400001.

Cheng J H, Hu L X, Wang T J, et al. Research progress of femtosecond laser-induced multifilament generation and regulation[J]. Chinese Journal of Lasers, 2023, 50(14): 1400001.

[7] Point G, Thouin E, Mysyrowicz A, et al. Energy deposition from focused terawatt laser pulses in air undergoing multifilamentation[J]. Optics Express, 2016, 24(6): 6271-6282.

[8] 刘伟伟, 薛嘉云, 苏强, 等. 超快激光成丝现象研究综述[J]. 中国激光, 2020, 47(5): 0500003.

Liu W W, Xue J Y, Su Q, et al. Research progress on ultrafast laser filamentation[J]. Chinese Journal of Lasers, 2020, 47(5): 0500003.

[9] Schroeder M C, Andral U, Wolf J P. Opto-mechanical expulsion of individual micro-particles by laser-induced shockwave in air[J]. AIP Advances, 2022, 12(9): 095119.

[10] Matthews M, Pomel F, Wender C, et al. Laser vaporization of cirrus-like ice particles with secondary ice multiplication[J]. Science Advances, 2016, 2(5): e1501912.

[11] Zeng Q W, Liu L, Ju J J, et al. Numerical investigation on the heat deposition characteristics of femtosecond laser pulses undergoing multiple filaments[J]. Physica Scripta, 2020, 95(8): 085605.

[12] 荀明娜, 尚滨鹏, 齐鹏飞, 等. 飞秒激光光丝空间特征的声学及荧光表征法: 对比研究[J]. 中国激光, 2023, 50(7): 0708008.

Xun M N, Shang B P, Qi P F, et al. Acoustic and fluorescence characterization of femtosecond laser filament spatial properties: comparative study[J]. Chinese Journal of Lasers, 2023, 50(7): 0708008.

[13] Rosenthal E W, Jhajj N, Larkin I, et al. Energy deposition of single femtosecond filaments in the atmosphere[J]. Optics Letters, 2016, 41(16): 3908-3911.

[14] Bychkov A S, Cherepetskaya E B, Karabutov A A, et al. Laser optoacoustic tomography for the study of femtosecond laser filaments in air[J]. Laser Physics Letters, 2016, 13(8): 085401.

[15] Wang L V, Hu S. Photoacoustic tomography: in vivo imaging from organelles to organs[J]. Science, 2012, 335(6075): 1458-1462.

[16] Potemkin F V, Mareev E I, Rumiantsev B V, et al. Two-dimensional photoacoustic imaging of femtosecond filament in water[J]. Laser Physics Letters, 2018, 15(7): 075403.

[17] 曾庆伟, 刘磊, 胡帅, 等. 基于多元线性阵列探测器的飞秒激光成丝光声图像重建[J]. 红外与激光工程, 2022, 51(8): 20210774.

Zeng Q W, Liu L, Hu S, et al. Femtosecond laser mercerization acoustic image reconstruction based on multivariate linear array detector[J]. Infrared and Laser Engineering, 2022, 51(8): 20210774.

[18] Rumiantsev B V, Mareev E I, Bychkov A S, et al. Three-dimensional hybrid optoacoustic imaging of the laser-induced plasma and deposited energy density under optical breakdown in water[J]. Applied Physics Letters, 2021, 118(1): 011109.

[19] 罗晓飞, 王波, 彭宽, 等. 基于聚焦声场模型的光声层析成像时间延迟快速校正反投影方法[J]. 物理学报, 2022, 71(7): 078102.

Luo X F, Wang B, Peng K, et al. Back-projection method with fast time-delay correction for photoacoustic tomography reconstruction based on a focused sound field model[J]. Acta Physica Sinica, 2022, 71(7): 078102.

[20] Xu M H, Wang L V. Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction[J]. Physical Review E, 2003, 67(5): 056605.

[21] Choi W, Oh D, Kim C. Practical photoacoustic tomography: realistic limitations and technical solutions[J]. Journal of Applied Physics, 2020, 127(23): 230903.

[22] Wahlstrand J K, Jhajj N, Rosenthal E W, et al. Direct imaging of the acoustic waves generated by femtosecond filaments in air[J]. Optics Letters, 2014, 39(5): 1290-1293.

[23] Mareev E I, Rumiantsev B V, Migal E A, et al. A comprehensive approach to the characterization of the deposited energy density during laser-matter interactions in liquids and solids[J]. Measurement Science and Technology, 2020, 31(8): 085204.

[24] Uryupina D S, Bychkov A S, Pushkarev D V, et al. Laser optoacoustic diagnostics of femtosecond filaments in air using wideband piezoelectric transducers[J]. Laser Physics Letters, 2016, 13(9): 095401.

[25] Bergé L, Skupin S, Lederer F, et al. Multiple filamentation of terawatt laser pulses in air[J]. Physical Review Letters, 2004, 92(22): 225002.

[26] Heritier J M. Electrostrictive limit and focusing effects in pulsed photoacoustic detection[J]. Optics Communications, 1983, 44(4): 267-272.

[27] Pramanik M. Improving tangential resolution with a modified delay-and-sum reconstruction algorithm in photoacoustic and thermoacoustic tomography[J]. Journal of the Optical Society of America A, 2014, 31(3): 621-627.

[28] Warbal P, Pramanik M, Saha R K. Impact of sensor apodization on the tangential resolution in photoacoustic tomography[J]. Journal of the Optical Society of America A, 2019, 36(2): 245-252.

[29] Pushkarev D, Mitina E, Shipilo D, et al. Transverse structure and energy deposition by a subTW femtosecond laser in air: from single filament to superfilament[J]. New Journal of Physics, 2019, 21(3): 033027.

[30] Sasoh A, Ohtani T, Mori K. Pressure effect in a shock-wave-plasma interaction induced by a focused laser pulse[J]. Physical Review Letters, 2006, 97(20): 205004.

[31] Liu W, Gravel J-F, Théberge F, et al. Background reservoir: its crucial role for long-distance propagation of femtosecond laser pulses in air[J]. Applied Physics B, 2005, 80(7): 857-860.

[32] 李骥, 张旻, Bogdan P. 空气耦合超声换能器声场的时域计算方法[J]. 无损检测, 2020, 42(5): 59-62.

Li J, Zhang M, Bogdan P. An approach to compute transient acoustic field radiated by air-coupled transducer[J]. Nondestructive Testing, 2020, 42(5): 59-62.

[33] 孙正, 孙慧峰. 生物光声层析成像中超声探测器特性对图像质量的影响及解决方法[J]. 中国生物医学工程学报, 2021, 40(6): 731-742.

Sun Z, Sun H F. Influence of ultrasonic detector characteristics on image quality in biological photoacoustic tomography and its solution[J]. Chinese Journal of Biomedical Engineering, 2021, 40(6): 731-742.

[34] Kalva S K, Pramanik M. Experimental validation of tangential resolution improvement in photoacoustic tomography using modified delay-and-sum reconstruction algorithm[J]. Journal of Biomedical Optics, 2016, 21(8): 086011.

曾庆伟, 刘磊, 胡帅, 李书磊, 赵世军. 光声层析重建飞秒光丝二维横向图像仿真研究[J]. 光学学报, 2024, 44(12): 1201011. Qingwei Zeng, Lei Liu, Shuai Hu, Shulei Li, Shijun Zhao. Numerical Study on Photoacoustic Tomography Reconstruction of Two-Dimensional Transverse Image of Femtosecond Laser Filaments[J]. Acta Optica Sinica, 2024, 44(12): 1201011.