高能激光热晕效应数值模拟方法的适用范围分析

热晕效应是影响高功率激光光束质量的重要因素之一，对这一效应进行合理的仿真，有利于高能激光的应用。针对目前模拟高能激光热晕效应的微扰法、积分法和相位屏法适用条件不清晰的问题，对三种数值模拟方法进行了系统的对比，并结合实验数据确定了每种方法最为合适的适用范围：广义畸变参数N<3时，选用积分法；3<N<4.8时，选用微扰法；N>4.8时，选用相位屏法。此外，利用可编程的液晶空间光调制器实现了热晕相位畸变的实验室模拟，得到的实验室模拟结果与数值仿真结果高度吻合。合理的数值模拟及实验室模拟对于强激光的实际应用具有重要的参考价值。

Abstract

Objective

Thermal blooming severely reduces beam quality and limits the efficiency of high-energy laser transmission in the atmosphere. Therefore, it is important to study the thermal blooming of high-energy laser propagation in the atmosphere systematically. Numerical simulation of thermal blooming is beneficial for the application of high-energy lasers. Numerical simulation methods for thermal blooming include perturbation, integral, and phase screen methods. However, most existing studies used only one of them for simulation, which led to obvious errors under some conditions. In this study, we compared the results of three methods and selected an appropriate numerical simulation range for each method based on the reported experimental results. In addition, field experiments for measuring the transmission of high-energy lasers in the atmosphere are complicated, and the experimental conditions are uncontrollable. Therefore, using a liquid crystal spatial light modulator (LC-SLM) based on the transmission characteristics of high-power lasers in the atmosphere is a valuable and cost-saving method for laboratory simulations of thermal blooming. We extracted the thermal blooming distortion phase based on the principle of the numerical simulation method and then used an LC-SLM to simulate thermal blooming in the laboratory. The experimental simulation results are consistent with the numerical simulation results.

Methods

In this study, numerical simulations of thermal blooming were performed using the perturbation, integral, and phase screen methods. First, the results of different numerical simulation methods for different transmission distances were compared. The results obtained using the different numerical simulation methods were significantly different, even under identical conditions. The relative peak intensity I_REL (ratio of maximum received light intensity to maximum initial light intensity) was measured as a function of the generalized distortion parameter N in the reported reference. A reasonably fitted curve was then selected as the reference data. The reasonable ranges for each method were determined by comparing the three numerical simulation results with the reference results. Subsequently, the thermal blooming distortion phase was extracted according to the applicable numerical simulation method under the conditions set for the laboratory experiments. Finally, laboratory experiments on thermal blooming were conducted using LC-SLM.

Results and Discussions

The relationship between the normalized peak intensity I_REL generated by the three numerical simulation methods and the generalized distortion parameter N is compared with the reference data at different transmission distances (Fig. 3). The errors of the perturbation and integral methods are small when N<4.8 (the errors of the two methods are approximate), and the error of the phase screen method is small when N>4.8. The relationship between the normalized peak intensity I_REL generated by the three numerical simulation methods and the generalized distortion parameter N is compared with the reference line at different initial laser powers (Fig. 4). The error of the perturbation method is small when N<3, that of the integral method is small when 3<N<4.8, and that of the phase screen method is small when N>4.8. The relationship between the normalized peak intensity I_REL generated by the three numerical simulation methods and the generalized distortion parameter N is compared with the reference line at different wind speeds (Fig. 5). The error of the perturbation method is small when N<2.8, that of the integral method is small when 2.8<N<4.8, and that of the phase screen method is small when N>4.8.

Conclusions

In this study, various ranges of generalized distortion parameter N applicable to each numerical simulation method are selected by comparing the error between the normalized peak intensity I_REL generated by the three numerical simulation methods and the reference line under different setting conditions. When the generalized distortion parameter N is less than 3, the error of the integral method is the smallest. When the generalized distortion parameter N ranges from 3 to 4.8, the error of the perturbation method is the smallest. When the generalized distortion parameter N is greater than 4.8, the error of the phase screen method is the smallest. Additionally, laboratory experiments are performed using LC-SLM. The phase of the thermal blooming distortion is accurately extracted and added to the phase modulator, and its effect of thermal blooming distortion is recorded using a CCD. The experimental results are in good agreement with the simulation results, confirming the feasibility of this experimental method. This work proposes a quantitative, accurate, programmable, and easily repeatable laboratory simulation device that provides an effective means for laboratory evaluations of the atmospheric transmission of high-energy lasers.

1　引言

高能激光在大气中传输时与大气相互作用产生的热晕效应严重限制了它在大气中的传输效率^［1-5］。热晕效应会使传输激光能量耗散以及光束质量下降，限制了高能激光的应用与发展^［6-13］。因此，研究高能激光在大气中传输的热晕效应具有重要意义^［14-17］。

高能激光的大气传输实验通常需要激光进行长距离传输，实验过程相当复杂，极为耗费人力物力，且实验条件不可控，进行室外实验的难度较大。因此，基于高功率激光在大气中的传输特性进行数值模拟以及实验室模拟具有重要的研究价值。近20年来，国内外研究人员对激光传输过程中的热晕效应进行了大量研究^［18-27］。2005年，Mahdieh等^［28］研究了环形CO₂激光束在高压CO₂气体中传输时的热晕效应，通过对比发现根据波动光学理论进行模拟得到的结果非常精确，特别是在小菲涅耳数情况下，而根据几何光学理论进行模拟得到的结果只对高菲涅耳数有效。2009年，姚友群^［29］利用微扰法模拟了有风情况下激光在大气中传输时的稳态热晕效应以及脉冲激光的瞬态热晕效应。2013年，Ji等^［30］研究了艾里光束在大气中传播的热晕效应，分析了艾里光束在大气中传输时光强分布、质心位置和均方光束宽度的变化。研究结果表明：热晕效应不仅会导致艾里光束中心波瓣的中心凹陷，还会导致中心波瓣光斑扩展和光强减小；与中心波瓣相比，旁瓣较少受到热晕效应的影响，这使得旁瓣的强度最大值可能大于中心波瓣的强度最大值。2016年，Shlenov等^［31］研究了CO₂激光在大气中传输时聚焦能力受风速的影响，结果显示：风速可以减弱热晕效应的影响，从而提高光束聚焦能力。2020年，吴书云^［32］采用数值模拟的方法对比了高斯光束、平顶光束、有中心遮拦的平顶光束在均匀大气传输过程中产生的稳态热晕效应，结果表明：在同样的广义热畸变参数N下，高斯光束的热晕效应最严重，平顶光束次之，空心光束的热晕效应最小。2021年，赵璐^［33］采用相位屏法研究了涡旋光束在大气中传输时传播距离、初始功率、横向风速、大气吸收系数、初始半径和波长对稳态热晕效应的影响。2022年，吴祥议政等^［34］研究了高功率激光在空气、氦气及氮气中传输时的热晕效应，发现氮气可以抑制热晕效应。

前述文献中所采用的研究热晕效应的数值模拟方法可以分为微扰法、积分法和相位屏法，这几种方法都各有优缺点以及一定的适用范围。在目前已有的研究中，人们大多只使用其中的一种方法进行各种条件下的模拟，没有考虑其适用范围的问题。鉴于此，笔者根据文献报道的实验数据，对这三种方法进行比较，尝试给出它们的适用范围。同时，笔者利用可编程的液晶空间光调制器（LC-SLM）实现热晕效应的实验室模拟，这种方法具有定量、可重复的优点。

本文第二节介绍稳态热晕效应的数值模拟，包括三种不同模拟方法的理论模型以及不同影响因素下的模拟结果；第三节将模拟结果与文献报道的实验结果进行对比分析，从而确定每种模拟方法的有效仿真区间；第四节介绍本次实验所用的液晶热晕模拟系统，并对实验结果进行讨论。最后，给出了本研究的结论。

2　热晕效应数值模拟

2.1　热晕效应数值模拟理论模型

微扰法将热晕效应产生的畸变当作一种简单的扰动项加载到原始未产生畸变的传输光束上，是一种几何光学方法，忽略了衍射效应的影响。若初始激光束为高斯光束，则激光传输稳态热晕效应的一阶微扰解^［25］可以表示为

\{\begin{array}{l} \frac{I (x, y, z)}{I_{u} (x, y, z) e x p (- α z)} = e x p \{- N_{c} \{2 (\frac{x}{a}) e x p (- \frac{x^{2} + y^{2}}{a^{2}}) + \frac{\sqrt[]{π}}{2} e x p (- \frac{y^{2}}{a^{2}}) (1 - 4 \frac{y^{2}}{a^{2}}) [1 + e r f (\frac{x}{a})]\}\}, \\ N_{c} = (- \frac{2 d n}{d T} \frac{I_{0} z}{n_{0} ρ c_{p} v a}) [1 - \frac{1 - e x p (- α z)}{α z}] \approx (- \frac{d n}{d T} \frac{I_{0} z^{2} α}{n_{0} ρ c_{p} v a}), α z ≪ 1, \end{array}

（1）

式中： $N_{c}$ 为热晕畸变参数； $d n / d T$ 为压强不变时大气折射率随温度的变化率； $I_{0}$ 为初始激光强度； $z$ 为传输距离； $n_{0}$ 为大气折射率； $ρ$ 为大气密度； $c_{p}$ 为比定压热容； $v$ 为风速； $α$ 为大气吸收系数； $a$ 为 $1 / e$ 光束半径；erf（·）为误差函数；I_u为初始激光光强分布。

积分法建立于惠更斯-菲涅耳原理之上。激光在各向同性、均匀、无源介质中传输的菲涅耳衍射积分方程^［35］为

u (x, y, z) = \frac{e x p (i k z)}{i λ z} e x p [i \frac{k}{2 z} (x^{2} + y^{2})] \iint \{u_{0} (x_{0}, y_{0}, 0) e x p [i \frac{k}{2 z} (x_{0}^{2} + y_{0}^{2})]\} e x p [- i \frac{k}{z} (x x_{0} + y y_{0})] d x_{0} d y_{0}

，（2）

式中： $λ$ 为激光波长；i为虚部单位；k为波数；x₀和y₀表示物平面的横纵坐标。积分法可以等效为：观察面上某一点的场强 $u (x, y, z)$ 与激光器发射初始波面上场强为 $u_{0} (x_{0}, y_{0}, 0)$ 且相位为 $\frac{k}{2 z} (x_{0}^{2} + y_{0}^{2})$ 的点光源的二维傅里叶变换成正比。在考虑热晕效应且有风的情况下，激光传输积分表达式^［36］为

\begin{array}{l} u (x, y, z) = \frac{e x p (i k z)}{i λ z} e x p [j \frac{k}{2 z} (x^{2} + y^{2})] \iint \{u_{0} (x_{0}, y_{0}, 0) e x p \{i [δ ψ + \frac{k}{2 z} (x_{0}^{2} + y_{0}^{2})]\}\} \cdot \\ e x p [- i \frac{k}{z} (x x_{0} + y y_{0})] d x_{0} d y_{0}, (3) \end{array}

其中 $δ ψ$ ^［36］可以表示为

δ ψ = - \frac{k (n_{0} - 1) (γ - 1) α z P}{2 \sqrt[]{π} γ p v a} e x p (- \frac{y_{0}^{2}}{a^{2}}) [1 + e r f (\frac{x_{0}}{a})]

，（4）

式中： $k$ 为波数； $γ = c_{p} / c_{V} \approx 1.4$ ， $c_{V}$ 为比定容热容； $P$ 为激光初始功率； $p$ 为大气压强。

相位屏法^［33］：在物与像之间插入多个像屏，对于相邻的两个屏，采用傅里叶变换法求出第二个屏上的光场分布，以此类推，得到最终像屏上的光场分布。对相邻两屏之间的光场分布进行求解的傅里叶变换法可以理解为：波传播半个步长后加上介质的作用传播下半个步长，求解真空波动方程，便可得到一个步长上的场。对光场进行傅里叶变换得到的光波表达式为

φ_{m n} (z, t) = \frac{1}{M^{2}} \sum_{p = 0}^{M - 1} \sum_{q = 0}^{M - 1} φ (p Δ x, q Δ y, z) e x p [\frac{- i 2 π (m p + n q)}{M}]

，（5）

式中：m和n分别为二维离散傅里叶变换中域矩阵的行数和列数；M为相位屏网格数。光波传播半个步长表示为

φ_{m n} (z + \frac{Δ z}{2}, t) = φ_{m n} (z, t) e x p [- \frac{i}{2 k} {(\frac{2 π}{L})}^{2} (m^{2} + n^{2}) \frac{Δ z}{2}]

，（6）

式中：L为相位屏宽度。将光波转化成光场，即

φ (p Δ x, q Δ y, z + \frac{Δ z}{2}, t) = \sum_{m = - \frac{M}{2}}^{\frac{M}{2} - 1} \sum_{n = - \frac{M}{2}}^{\frac{M}{2} - 1} φ_{m n} (z + \frac{Δ z}{2}, t) e x p [\frac{- i 2 π (m p + n q)}{M}]

，（7）

加入介质作用后，光场表示为

φ' (p Δ x, q Δ y, z + \frac{Δ z}{2}, t) = φ (p Δ x, q Δ y, z + \frac{Δ z}{2}, t) e x p (\frac{i k}{2} δ \bar{ε} Δ z)

，（8）

其中 $δ \bar{ε} = \frac{1}{Δ z} \int_{z^{n}}^{z^{n} + Δ z} δ ε (z') d z'$ ， $δ ε = 0.46 ρ_{1}$ ， $ρ_{1} = \frac{(γ - 1) α}{c_{s}^{2} v} \int_{- \infty}^{x} I (x', y, z) d x'$ ， $c_{s}^{2}$ 为声速的平方。将光场表达式转换成光波，即

φ'_{m n} (z, t) = \frac{1}{M^{2}} \sum_{p = 0}^{M - 1} \sum_{q = 0}^{M - 1} φ' (p Δ x, q Δ y, z + \frac{Δ z}{2}, t) e x p [\frac{- i 2 π (m p + n q)}{M}]

。（9）

传播下半个步长可表示为

φ'_{m n} (z + Δ z, t) = φ'_{m n} (z, t) e x p [- \frac{i}{2 k} {(\frac{2 π}{L})}^{2} (m^{2} + n^{2}) \frac{Δ z}{2}]

，（10）

从而可得出传播一个步长的光场分布为

φ' (p Δ x, q Δ y, z + Δ z, t) = \sum_{m = - \frac{M}{2}}^{\frac{M}{2} - 1} \sum_{n = - \frac{M}{2}}^{\frac{M}{2} - 1} φ'_{m n} (z + Δ z, t) e x p [\frac{- i 2 π (m p + n q)}{M}]

。（11）

参考以上一个步长的求解过程进行类推，即可得出整个传输过程。

2.2　热晕效应数值模拟结果比较

存在横向风时，激光传输热晕效应可以看作是由强制对流产生的，即空气沿风向穿过激光束时，空气中的分子不断吸收激光能量，这部分空气的温度越来越高，导致折射率降低，光束偏折向折射率大的一方，从而形成月牙形热晕畸变。设波长λ=10.6 μm，激光发射功率P=50 kW，光束半径a=0.2 m，大气吸收系数α=6.5×10^-5 m^-1，风速v=5 m/s，传输距离z=5 km，发射激光光束为近准直高斯光束。下面将讨论不同模拟方法下改变各参数对畸变产生的影响。

在上述设定条件下，分别取传输距离为1、3、5、10 km，利用微扰法、积分法、相位屏法分别模拟归一化光强（接收光强/接收光强最大值）等值线分布，模拟结果如图1所示。随着传输距离不断增加，畸变参数N_c不断增大，热晕效应的影响增强，光束逐渐变为月牙形。在条件相同的情况下，利用不同的方法进行模拟的结果各不相同：当传输距离较远时，微扰法的模拟结果没有明显的扩束现象，积分法的模拟结果与实际吻合度不高；当传输距离较近时，相位屏法的模拟结果展现出来的畸变效果相对较强。类似地，在不同初始功率和风速条件下，三种方法的仿真结果也有明显的区别。

图 1. 三种方法模拟的不同传输距离下的光强等值线图。（a1）（b1）（c1）（d1）微扰法模拟结果；（a2）（b2）（c2）（d2）积分法模拟结果；（a3）（b3）（c3）（d3）相位屏法模拟结果

Fig. 1. Contour diagrams of light intensity simulated by the three methods at different transmission distances. (a1)(b1)(c1)(d1) Simulation results of perturbation method; (a2)(b2)(c2)(d2) simulation results of integral method; (a3)(b3)(c3)(d3) simulation results of phase screen method

下载图片查看所有图片

综上可知：随着传输距离、初始激光功率以及吸收系数增大，热晕效应增强，激光传输效率降低；随着风速、光束半径以及激光波长增大，热晕效应减弱，激光传输效率提高。另外可以看出，在相同的参数下，不同数值模拟方法的模拟结果具有较大差异。接下来，本文将讨论不同数值模拟方法的适用范围。

3　三种模拟方法的适用范围

在以上三种模拟方法中，微扰法最为简单，但当畸变参数达到一定值时微扰法将会不准确，同时微扰法无法展现光斑扩散现象；积分法的准确度相比微扰法得到了一定的提升，但当畸变参数值变大时，尤其是当传输距离增大从而导致畸变参数增大时，模拟结果也会不准确；相位屏法，目前广泛使用的一种方法，在畸变值很小时，呈现出来的畸变效果十分明显。故而，为每种方法找到其最为合适的使用范围，将会提高数值模拟的准确性。图2是Gebhardt^［26］于1990年发表的“Twenty-five years of thermal blooming: an overview”论文中展示的实验结果与数值模拟结果中的I_REL与N的关系对比图，其中I_REL指的是归一化峰值强度（接收光强最大值/初始光强最大值），N指的是广义畸变参数。这张图片显示实验结果与数值模拟结果十分吻合，误差在0.3左右。因实验条件受限，笔者参考图2中的实线来确定每种模拟方法的适用范围。该实线的表达式^［26］为

I_{R E L} = \frac{1}{1 + 0.0625 N^{2}}

。（12）

图 2. 实验和数值模拟的I_REL与N的关系图^［26］

Fig. 2. Relative peak intensity I_REL dependence on the distortion parameter N for experimental and simulated results^[26]

下载图片查看所有图片

为方便与参考线进行对照，选用广义畸变参数N为横坐标。下文将绘制不同数值模拟方法下归一化峰值强度I_REL与N的关系图。广义畸变参数的表达式^［26］为

N = N_{c} f (N_{E}) q (N_{F}) s (N_{ω})

，（13）

其中，

\{\begin{array}{l} q (N_{F}) = [\frac{2 N_{F}^{2}}{N_{F} - 1}] [1 - \frac{l n N_{F}}{N_{F} - 1}] \\ f (N_{E}) = \frac{2}{N_{E}^{2}} [N_{E} - 1 + e x p (- N_{E})] \\ s (N_{ω}) = \frac{2}{N_{ω}^{2}} [(N_{ω} + 1) l n (N_{ω} + 1) - N_{ω}] \end{array}

，（14）

式中： $f (N_{E}) 、 q (N_{F}) 、 s (N_{ω})$ 分别为消光数修正因子、菲涅耳数修正因子、光束扫描修正因子；N_c为准直光束稳态热晕畸变参数。N_c以及N_E、N_F、N_ω的计算公式分别为

N_{c} = \frac{16 \sqrt[]{2} (- n_{T}) α P L^{2}}{π n_{0} ρ c_{p} v D^{3}},

（15）

N_{E} = (α + s) z,

（16）

N_{F} = \frac{k a^{2}}{z},

（17）

N_{ω} = \frac{ω z}{v_{0}},

（18）

式中：n_T表示压强不变时大气折射率随温度的变化率；D为激光光束直径。

为了更直观地观察每种方法的适用范围，本文引入误差这一参考量，其可表示为

\{\begin{matrix} I_{e r r o r 1} = I_{p e r t u r b a t i o n} - I_{E m p i r i c a l} \\ I_{e r r o r 2} = I_{p h a s e s c r e e n} - I_{E m p i r i c a l} \end{matrix}

，（19）

式中： $I_{e r r o r}$ 表示模拟结果与实验结果的误差； $I_{p e r t u r b a t i o n}$ 表示使用微扰法模拟的归一化峰值强度； $I_{p h a s e s c r e e n}$ 表示使用相位屏法模拟的归一化峰值强度； $I_{E m p i r i c a l}$ 表示参考线的归一化峰值强度。

3.1　传输距离的影响

取传输距离为1~50 km，步长为0.1 km，其他条件保持不变。图3（a）所示为不同传输距离下三种数值模拟方法模拟的热晕效应的归一化峰值强度（I_REL）与广义畸变参数之间的关系，图3（b）所示为不同传输距离下三种数值模拟方法模拟的热晕效应的归一化峰值强度与参考线的误差（I_error）分析图。将图3（a）与图3（b）结合起来观察可以发现：当N<4.8时，选用微扰法或者积分法（两种方法的误差比较接近）比较合适；当N>4.8时，选用相位屏法比较合适。

图 3. 不同传输距离下三种数值模拟方法的适用范围。（a）I_REL与N的关系；（b）I_error与N的关系

Fig. 3. Application range of three numerical simulation methods under different transmission distances. (a) I_REL as a function of N; (b) I_error as a function of N

下载图片查看所有图片

3.2　初始激光功率的影响

取初始激光功率为1~300 kW，步长为0.5 kW，其他条件保持不变。图4（a）所示为不同初始激光功率下三种数值模拟方法模拟的热晕效应的归一化峰值强度与广义畸变参数的关系，图4（b）所示为不同初始激光功率下三种数值模拟方法模拟的热晕效应的归一化峰值强度与参考线的误差分析图。将图4（a）与图4（b）结合起来观察可以发现：当N<3时，选用积分法比较合适；当3<N<4.8时，选用微扰法比较合适；当N>4.8时，选用相位屏法比较合适。从图4（b）中可以观察到N>25时积分法模拟的误差更小，但最终没有在这个区间选择该方法。这是由于初始激光功率增大到一定值时，尽管积分法模拟的热晕效应的峰值强度较高，但它的形状不再呈现月牙形，即此时的畸变参数已经超出了积分法可模拟的范围。

图 4. 不同初始激光功率下三种数值模拟方法的适用范围选取。（a）I_REL与N的关系；（b）I_error与N的关系

Fig. 4. Application range of three numerical simulation methods under different initial laser powers. (a) I_REL as a function of N; (b) I_error as a function of N

下载图片查看所有图片

3.3　风速的影响

取风速为0.5~30 m/s，步长为0.2 m/s，其他条件保持不变。图5（a）所示为不同风速下三种数值模拟方法模拟的热晕效应的归一化峰值强度与广义畸变参数的关系，图5（b）为不同风速下三种数值模拟方法模拟的热晕效应的归一化峰值强度与参考线的误差分析图。将图5（a）与图5（b）结合起来观察可以发现：当N<2.8时，选用积分法比较合适；当2.8<N<4.8时，选用微扰法比较合适；当N>4.8时，选用相位屏法比较合适。

图 5. 不同风速下三种数值模拟方法的适用范围。（a）I_REL与N的关系；（b）I_error与N的关系

Fig. 5. Application range of three numerical simulation methods under different wind speeds. (a) I_REL as a function of N; (b) I_error as a function of N

下载图片查看所有图片

尽管每种影响因素对应的N值选用的区间不同，但数值上都是比较接近的，考虑到在不同模型模拟曲线相交之处所对应的N值附近选取N值区间时三种模型的误差较小，得出以下结论：当广义畸变参数N<3时，选用积分法；当3<N<4.8时，选用微扰法；当N>4.8时，选用相位屏法。这样的区间选取法可使数值模拟结果与参考线相比存在0.2以内的误差，与实验结果的误差约在0.1以内，证明了该方法的可行性。经分析，此处模拟误差来源可能包括以下三个方面：1）模拟过程中设定大气吸收系数α为固定值，忽略了激光在大气中传输时引起的大气温度变化会导致大气密度改变，从而引起大气吸收系数变化；2）实验数据无规律地分布在所选用的参考线附近，实验数据与参考线之间本身就存在一定的误差；3）激光在大气中传输的实验环境过于复杂，存在不可控因素，无法像数值模拟一样给出确定的每个量。

4　实验结果与分析

由于激光传输实验所耗人力物力巨大，因此笔者基于实验室已有的器材，选用LC-SLM进行实验室模拟。搭建的实验光路如图6所示。激光器为氦氖内腔激光器，波长为632.8 nm。LC-SLM的光束输出类型为反射式，分辨率为512 pixel×512 pixel，灰度阶数为256阶8位。透镜1的焦距为20 mm，透镜2的焦距为400 mm。

图 6. 基于LC-SLM的热晕效应实验光路图

Fig. 6. Optical path of thermal blooming effect experiment based on LC-SLM

下载图片查看所有图片

选用液晶空间光调制器（LC-SLM）进行热晕效应实验，实验过程中最主要的是进行热晕畸变相位的提取。从理论上来说，热晕效应数值模拟的三种方法（微扰法、积分法、相位屏法）均可进行热晕畸变相位的提取。下面将分别介绍三种数值模拟方法进行热晕畸变相位提取的过程。

基于微扰法的数值模拟方法的相位提取可以参考文献［26］，该文提到热晕相位可以表示为

表 1. 光斑参数的理论和实验值对比

Table 1. Theoretical and experimental values comparison of light spot parameters

z /km	m_theoretical	m_experimental	x_theoretical	x_experimental
1	242	32	256.0611	258.6361
3	226	30	252.1005	256.9763
5	205	28	244.3762	254.9481
10	159	22	225.8691	251.6324

查看所有表

φ_{B G} (x, y, z) = - \frac{Δ φ_{G}}{2} e x p [- {(\frac{y}{a})}^{2}] [1 + e r f (\frac{x}{a})]

，（20）

其中，

Δ φ_{G} = \frac{1}{2 \sqrt[]{π}} N_{D}

，（21）

N_{D} = \frac{4 \sqrt[]{2} (- n_{T}) k α P z}{ρ c_{p} v D}

。（22）

基于积分法的数值模拟方法的相位提取可以参考本文第2节介绍的积分法的基本原理。对于热晕效应来说，激光传输过程中介质吸收热量后折射率发生的变化为 $δ n$ ，该扰动的最终影响可认为是在激光束初始波面上加一个相位扰动项 $δ ψ$ ，使光束传播发生变化。那么，此时的相位扰动项 $δ ψ$ 即为积分法中提取出的畸变相位，可以表示为

\begin{array}{l} δ ψ = - \frac{k (n_{0} - 1) (γ - 1) α z P}{2 \sqrt[]{π} γ p v a} e x p (- \frac{y_{0}^{2}}{a^{2}}) \cdot \\ [1 + e r f (\frac{x_{0}}{a})] \end{array}

，（23）

式中： $γ = c_{p} / c_{V} \approx 1.4$ 。

基于相位屏法的数值模拟方法的相位提取可以参考文献［37］，相邻相位屏之间的相位畸变可以表示为

Δ φ = k δ \bar{ε} Δ z / 2

，（24）

其中， $δ \bar{ε} = \frac{1}{Δ z} \int_{z^{n}}^{z^{n} + Δ z} δ ε (z') d z'$ ， $δ ε = 0.46 ρ_{1}$ ， $ρ_{1} = \frac{(γ - 1) α}{c_{s}^{2} v} \int_{- \infty}^{x} I (x', y, z) d x'$ 。整个传播过程中的相位畸变为每个步长下相位畸变累加的结果。

设波长λ=10.6 μm，激光发射功率P=50 kW，光束半径a=0.25 m，大气吸收系数α=6.5×10^-5 m^-1，风速v=5 m/s，发射激光光束为高斯光束。传输距离为1、3、5、10 km时，N值分别为0.7047、2.0192、3.1581、5.6045。根据第3节得到的结论，当z为1 km和3 km时使用积分法进行实验，当z为5 km时使用微扰法进行实验，当z为10 km时使用相位屏法进行实验。实验结果如图7所示。从图中可以明显观察到：随着传输距离增大，热晕畸变增强，光束由最初的圆斑逐渐变为月牙形。

图 7. 热晕实验结果

Fig. 7. Thermal blooming experimental results

下载图片查看所有图片

在实验过程中，笔者记录了光束横向长度及光束质心位置，并将其与数值模拟的结果进行对比，如表1所示。其中，m_theoretical和x_theoretical分别表示数值模拟结果中光束的横向长度以及质心位置的横坐标，m_experimental和x_experimental分别表示实验室模拟结果中光束的横向长度以及质心位置的横坐标，须注意的是长度均由像素点的个数来表示。

使用Pearson积矩相关系数^［38］对表1中的数据进行计算分析，计算公式为

ρ_{x, y} = \frac{c o v (x, y)}{σ_{x} σ_{y}}

，（25）

式中：cov（·）是变量之间的协方差；σ是变量的标准差。将m_theoretical和m_experimental代入式（25）得到相关系数为0.9858，将x_theoretical和x_experimental代入式（25）得到相关系数为0.9896，相关系数的绝对值很接近1，表明实验室模拟结果与数值模拟结果极为接近。

5　结论

本文着重介绍了研究热晕效应常用的三种数值模拟方法，分析了这三种方法的优缺点，并参考实验结果，以广义畸变参数N为参考对象为每种方法找到了其最为合适的适用范围，即：当N<3时，选用积分法；当3<N<4.8时，选用微扰法；当N>4.8时，选用相位屏法。这样得到的模拟结果相对于参考线的误差在0.2之内，而相较于实验结果的误差会更小。除此之外，还进行了实验室模拟实验。选用LC-SLM进行模拟实验，将热晕畸变的相位添加到相位调制器上，利用相机观察热晕畸变的效果。结果显示，实验室模拟结果与数值模拟结果的吻合度较高，验证了实验室模拟方法的可行性。

本工作提供了一种定量、准确、可编程、易重复的实验室模拟装置，为高功率激光的大气传输实验室评价提供了有效手段。

参考文献

[1] Van Zandt N R, Fiorino S T, Keefer K J. Enhanced, fast-running scaling law model of thermal blooming and turbulence effects on high energy laser propagation[J]. Optics Express, 2013, 21(12): 14789-14798.

[2] Banakh V A, Falits A V. Efficiency of combined beam focusing under thermal blooming[J]. Atmospheric and Oceanic Optics, 2014, 27(3): 211-217.

[3] Sprangle P, Hafizi B, Ting A, et al. High-power lasers for directed-energy applications[J]. Applied Optics, 2015, 54(31): F201-F209.

[4] 李晓庆, 季小玲. 大气湍流和热晕效应对列阵合成激光传输特性和光束质量影响的理论研究进展[J]. 强激光与粒子束, 2023, 35(4): 041007.

Li X Q, Ji X L. Theoretical research progress on the influence of atmospheric turbulence and thermal blooming on characteristics and beam quality of laser array beams propagating in the atmosphere[J]. High Power Laser and Particle Beams, 2023, 35(4): 041007.

[5] 徐梦敏, 李晓庆, 唐荣, 等. 风控热晕对双模涡旋光束大气传输的轨道角动量和相位奇异性的影响[J]. 物理学报, 2023, 72(16): 164202.

Xu M M, Li X Q, Tang R, et al. Influence of wind-dominated thermal blooming on orbital angular momentum and phase singularity of dual-mode vortex beams[J]. Acta Physica Sinica, 2023, 72(16): 164202.

[6] Li X Q, Ji X L. Effect of thermal blooming on the beam quality of truncated laser beams[J]. Optik, 2016, 127(20): 8350-8356.

[7] Wang L, Lin W B. The effect of wind directions on the thermal blooming of a laser beam propagating in the air[J]. Journal of Russian Laser Research, 2016, 37(4): 407-410.

[8] 邓凌, 李晓庆, 穆轶, 等. 风控热晕下椭圆激光光束质量的优化研究[J]. 中国激光, 2022, 49(4): 0405001.

Deng L, Li X Q, Mu Y, et al. Beam quality optimization of elliptical laser under wind-dominated thermal blooming[J]. Chinese Journal of Lasers, 2022, 49(4): 0405001.

[9] Wang L, Lin W B. The effect of the wind speed on the thermal blooming of a laser beam propagating in the rain[J]. Journal of Russian Laser Research, 2017, 38(5): 455-458.

[10] Yu L, Hu B B, Zhang Y X. Intensity of vortex modes carried by Lommel beam in weak-to-strong non-Kolmogorov turbulence[J]. Optics Express, 2017, 25(16): 19538-19547.

[11] Akers B F, Reeger J A. Numerical simulation of thermal blooming with laser-induced convection[J]. Journal of Electromagnetic Waves and Applications, 2019, 33(1): 96-106.

[12] Kanev F Y, Makenova N A, Lukin V P, et al. Adaptive correction of thermal distortions of multichannel laser radiation[J]. Atmospheric and Oceanic Optics, 2018, 31(3): 238-242.

[13] Ding Z L, Li X Q, Cao J Y, et al. Thermal blooming effect of Hermite-Gaussian beams propagating through the atmosphere[J]. Journal of the Optical Society of America A, 2019, 36(7): 1152-1160.

[14] Wu S Y, Luo X, Li X Y. Adaptive optics for reduction of thermal blooming effects by the phase compensation[J]. Journal of Russian Laser Research, 2020, 41(4): 413-423.

[15] Ahn K, Lee S H, Park I K, et al. Numerical simulation of high-energy laser propagation through the atmosphere and phase correction based on adaptive optics[J]. Journal of the Korean Physical Society, 2021, 79(10): 918-929.

[16] Lu L, Wang Z Q, Zhang P F, et al. Thermal blooming induced phase change and its compensation of a Gaussian beam propagation in an absorbing medium[J]. Optics Letters, 2021, 46(17): 4304-4307.

[17] Qiu D, Tian B Y, Ting H, et al. Mitigation of thermal blooming by rotating laser beams in the atmosphere[J]. Applied Optics, 2021, 60(27): 8458-8465.

[18] McLean E A, Sica L, Glass A J. Interferometric observation of absorption induced index change associated with thermal blooming[J]. Applied Physics Letters, 1968, 13(11): 369-371.

[19] Gebhardt F G, Smith D C. Effects of wind on thermal defocusing of CO2 laser radiation[J]. Applied Physics Letters, 1969, 14(2): 52-54.

[20] Smith D. Thermal defocusing of CO2 laser radiation in gases[J]. IEEE Journal of Quantum Electronics, 1969, 5(12): 600-607.

[21] Chan C H. Effective absorption for thermal blooming due to aerosols[J]. Applied Physics Letters, 1975, 26(11): 628-630.

[22] Fleck J A, Jr, Morris J R, Feit M D. Time-dependent propagation of high energy laser beams through the atmosphere[J]. Applied Physics, 1976, 10(2): 129-160.

[23] Gebhardt F G. High power laser propagation[J]. Applied Optics, 1976, 15(6): 1479-1493.

[24] Primmerman C A, Fouche D G. Thermal-blooming compensation: experimental observations using a deformable-mirror system[J]. Applied Optics, 1976, 15(4): 990-995.

[25] Smith D C. High-power laser propagation: thermal blooming[J]. Proceedings of the IEEE, 1977, 65(12): 1679-1714.

[26] Gebhardt F G. Twenty-five years of thermal blooming: an overview[J]. Proceedings of SPIE, 1990, 1221: 2-25.

[27] Vorob'ev V V. Thermal blooming of laser beams in the atmosphere[J]. Progress in Quantum Electronics, 1991, 15(1/2): 1-152.

[28] Mahdieh M H, Lotfi B. Two-dimensional-simulation of thermal blooming effects in ring pattern laser beams[J]. Optical Engineering, 2005, 44(9): 096001.

[29] 姚友群. 强激光在大气中传输的热晕效应研究[D]. 西安: 西安电子科技大学, 2009.

YaoY Q. Study on thermal halo effect of intense laser propagation in the atmosphere[D]. Xi'an: Xidian University, 2009.

[30] Ji X L, Eyyuboğlu H T, Ji G M, et al. Propagation of an Airy beam through the atmosphere[J]. Optics Express, 2013, 21(2): 2154-2164.

[31] Shlenov S A, Vasiltsov V V, Kandidov V P. Energy parameters of CO2 laser radiation focused in a turbulent atmosphere under wind-dominated thermal blooming[J]. Atmospheric and Oceanic Optics, 2016, 29(4): 324-330.

[32] 吴书云. 激光大气传输热晕效应及自适应光学校正的数值仿真和实验研究[D]. 成都: 中国科学院光电技术研究所, 2020.

WuS Y. Numerical simulation and experimental study on thermal blooming effect and adaptive optical correction of laser propagation in atmosphere[D]. Chengdu: Institute of Optics and Electronics, Chinese Academy of Sciences, 2020.

[33] 赵璐. 涡旋光束在大气中传输的热晕效应[D]. 太原: 太原科技大学, 2021. 10.1016/j.optlastec.2021.106982

ZhaoL. Thermal halo effect of vortex beam propagation in atmosphere[D]. Taiyuan: Taiyuan University of Science and Technology, 2021.

[34] Wu X Y Z, Xu J, Gong K L, et al. Theoretical and experimental studies on high-power laser-induced thermal blooming effect in chamber with different gases[J]. Chinese Physics B, 2022, 31(8): 086105.

[35] 宫鹏, 徐桂成, 张斯特, 等. 菲涅耳衍射及有限时域差分在光栅仿真中的应用研究[J]. 长春理工大学学报(自然科学版), 2022, 45(3): 1-5.

Gong P, Xu G C, Zhang S T, et al. Research of Fresnel diffraction and finite difference time domain in grating simulation[J]. Journal of Changchun University of Science and Technology (Natural Science Edition), 2022, 45(3): 1-5.

[36] 蔡邦维, 黄文龙, 吕百达, 等. CW稳态热晕的数值模拟[J]. 光学学报, 1996, 16(7): 948-951.

Cai B W, Huang W L, Lü B D, et al. Numerical simulation of CW steady-state thermal blooming with wind[J]. Acta Optica Sinica, 1996, 16(7): 948-951.

[37] Zhang Q, Hu Q L, Wang H Y, et al. High-precision calculation and experiments on the thermal blooming of high-energy lasers[J]. Optics Express, 2023, 31(16): 25900-25914.

[38] 盛骤, 谢式千, 潘承毅. 概率论与数理统计[M]. 4版. 北京: 高等教育出版社, 2008.

ShengZ, XieS Q, PanC Y. Probability and mathematical statistics[M]. 4th ed. Beijing: Higher Education Press, 2008.

张琪, 胡启立, 王红燕, 胡鸣, 徐星宇, 吴晶晶, 胡立发. 高能激光热晕效应数值模拟方法的适用范围分析[J]. 中国激光, 2024, 51(8): 0805003. Qi Zhang, Qili Hu, Hongyan Wang, Ming Hu, Xingyu Xu, Jingjing Wu, Lifa Hu. Analysis of Applicable Scope of Numerical Simulation Method for Thermal Blooming of High-Energy Laser[J]. Chinese Journal of Lasers, 2024, 51(8): 0805003.