0 引言

当激光束在大气中传输时，除受大气微粒的吸收、散射外还将受到大气湍流的影响，造成额外的光束扩展、光束瞬时中心的随机漂移以及光强的随机扰动^[1-4].由于大气湍流的随机特性，使我们不得不从统计的角度描述大气湍流的特征，而围绕大气湍流所建立的各种解析理论都是基于一定的假设或近似，尚未有一种统一的理论能够完整地描述各种大气条件下的湍流特性.

湍流的数值模拟方法，特别是以多相位屏模型为代表的分步传输(split-step)方法，是突破各种假设和近似，从而获得各种条件下大气湍流特征的有效方法，但因其模糊的数学表达式以及庞大的计算成本，使得其早期仅用于验证解析理论的准确性^[5-7].但随着近年来计算机技术的飞速发展，数值模拟方法已经摆脱计算速度慢的局面，使其成为了一种研究激光大气湍流传输问题的有效手段，因此数值模拟的精度自然成为了最重要的指标.数值模拟精度主要受相位屏生成精度的影响.针对相位屏的生成，国内外学者们已提出了很多不同的方法，其中具有代表性的有基于快速傅立叶变换(Fast Fourier Transform，FFT)的蒙特卡罗方法^[8]、Zernike多项式法^[9]等，其中基于FFT的相位屏生成方法因其快速的计算优势，成为了湍流相位屏模拟的首选方法.但是FFT方法具有一个严重的缺陷，即低频部分功率的缺失，使得相位屏在表征大尺度湍涡的特征上存在较大误差，因此要使用FFT方法生成的相位屏进行湍流传输模拟，必须对低频成分进行补偿.作者曾在文献[1]中提出一种基于优化的相位屏生成方法(Optimization-Based Mthod，OBM)^[1]，该模型首先利用FFT方法生成高频相位屏，然后对折射率谱的低频区域进行稀疏重采样，并借助优化算法找出使得相位屏模拟误差最小的采样点分布，通过这些采样点生成低频补偿屏，从而使得相位屏低频误差得到补偿.

文献[1]中，作者所提出的采样设置针对的是Von-Karman谱，这是推导解析理论时常用的谱模型，在数学上易于处理并且在极限状态下能够简化为Kolmogorov谱和Tatarskii谱^[10-11].然而这些谱模型仅在惯性区与真实的湍流折射率谱一致，而超出惯性区的部分是为了数学处理的方便而并非物理实际，实验表明，湍流折射率谱在接近1/l₀的高波数区存在一个小的“凸起”^[12]，使得谱的下降速率小于上述谱所预测的κ^-11/3(κ为波数).文献[13]指出折射率谱的“凸起”将会同样的反映在光束的各种统计量上，例如结构函数以及光束漂移等.因此要想实现对大气湍流的准确模拟，有必要将此“凸起”考虑进去，ANDREWS L C等提出了一种修正大气谱，通过分析准确的Hill数值谱，在传统Von-Karman谱的基础上进行了修正，使其既具有Hill谱所描述的“凸起”又具有类似于Von-Karman谱的解析表达式^[14].

本文试图通过之前提出的基于优化的相位屏模型实现具有修正大气谱特征的相位屏的高精度生成以及高斯光束的湍流传输模拟.首先根据修正大气谱的特征，改进了原始低频区域的采样设置，实现了基于修正大气谱的高精度相位屏的生成；然后，以生成的高精度相位屏为基础，采用split-step方法对高斯光束的湍流传输进行了数值模拟，对光束扩展和光束漂移进行了统计，最后与ANDREWS L C等利用该修正谱推导的解析结果^[14-15]进行了比较.

1 基于修正大气谱的相位屏模型

作者针对传统的Von-Karman谱，在由高频网格第一采样点所围成的正方形区域内选取了16个采样点，通过两个位置参数(f₁和f₂)改变采样点的分布，如图 1所示.

图 1. 低频采样点分布

Fig. 1. The distribution of low frequency samples

下载图片 查看所有图片

通过优化算法使得最终生成的相位屏在低频部分的最大相对误差降低至0.1%以下，但是文中采用的Von-Karman谱仅考虑了外尺度影响，这对于光束扩展以及光束漂移等低阶统计量的影响并不明显^[16]，但是对于闪烁指数这类高阶统计量来说，内尺度的存在极大地影响了闪烁指数的估计值，特别是在强湍流情况下^[17-18]，此外正如上文所提到的，折射率谱在高波数处的“凸起”也会影响到各类统计量的估计.因此单纯使用Von-Karman谱无法对激光束通过大气湍流后的全部特征进行模拟，必须选用更贴近实际的修正大气谱，如式(1)所示^[12]

式中，$C_{\rm n}^2$为大气折射率结构常数，κ为空间波数，κ_l=3.3/l₀，l₀为内尺度，κ₀=2π/L₀，L₀为外尺度，a₁=1.802，a₂=0.254，以Kolmogorov谱为尺度，对Tatarskii谱、Von-Karman谱以及修正大气谱进行归一化，结果如图 2所示.

图 2. 归一化折射率谱模型

Fig. 2. Scaled spectral models of refractive-index fluctuations

下载图片 查看所有图片

从图 2中可以看出，修正大气谱相对于其他两个谱来说，具有明显的“凸起”，此外对比Von-Karman谱与另外两个谱可以看出在更高的波数区，Tatarskii谱及修正大气谱由于内尺度的影响将会以更高的速率衰减.修正大气谱的“凸起”使得该谱在凸起前的低频部分具有比Kolmogorov谱更缓的衰减，在凸起后的高频部分具有更快的衰减，因此对修正大气谱进行采样时，不再像Von-Karman谱一样，绝大部分能量集中在第一采样点所包含的范围内，而是会过渡到第二采样点甚至更高频率的采样点所包含的区域.基于此，我们首先将低频采样区域由原先的第一采样点包围的区域扩大至第二采样点包围的区域，然后对第一、第二采样点所包围的区域进行独立的采样，每个区域的采样点设置与文献[1]中一致，总共32个采样点以及4个位置参数，分别为κ₁, κ₂∈[0, 2π/L]以及κ₃, κ₄∈[0, 4π/L]，采样点以空间波数表示，如图 3所示.

图 3. 修正大气谱采样点分布

Fig. 3. The samples distribution of modified atmosphere spectrum

下载图片 查看所有图片

通过对第一、第二采样点所包围的区域进行独立采样，可分别生成两个低频补偿屏，根据稀疏谱模型，这两个低频补偿屏可表示为

式中，

分别为第一采样点包围区域和第二采样点包围区域的相位谱，${r_0} = {(0.432{k^2}C_n^2Z)^{ - 3/5}}$为Fried参数，表示湍流的强弱，k=2π/λ为入射光波数，Z为传输距离，κ_xjk、κ_yjk和κ_xuv、κ_yuv分别为图 3(a)和图 3(b)所示采样点的空间波数坐标，Δκ_xjk、Δκ_yjk和Δκ_xuv、Δκ_yuv为各采样点的波数间隔，L_x和L_y为相位屏在x、y方向的尺寸，N_x和N_y分别为相位屏在x、y方向的采样点数，在此假设相位屏为方形屏，即L_x=L_y=L且N_x=N_y=N，m、n和m^′、n^′为整数，h(j, k)和h^′(j, k)为零均值、单位方差的复高斯随机矩阵.

将两个低频补偿屏与高频屏相加即可得到补偿后的精确相位屏，高频相位屏的生成可参见文献[1]，只需将文献[1]中式(3)中的Von-Karman谱变为修正大气谱即可，在此不再赘述.由上述分析可知，该补偿方案中存在4个自由参数，用于确定两个低频补偿屏的采样点分布，其寻优方法与文献[1]中一致，依然采用引力搜索算法，其中待优化参数为c₁, c₂, c₃, c₄∈[0, 100]，与位置参数的关系为${\kappa _1} = {c_1}\cdot\frac{{2{\rm \pi }}}{{100\cdotL}},{\kappa _2} = {c_2}\cdot\frac{{2{\rm \pi }}}{{100\cdotL}},{\kappa _3} = {c_3}\cdot\frac{{4{\rm \pi }}}{{100\cdotL}},{\kappa _4} = {c_4}\cdot\frac{{4{\rm \pi }}}{{100\cdotL}}$，样本维度由2维变为4维，理论结构函数变为修正大气谱的理论结构函数，其他参数保持不变.根据Andrews的推导，修正大气谱的结构函数可表示为^[19]

式中，r为空间中两点间的距离.

2 相位屏生成及高斯光束传输模拟

相位屏的固有参数r₀，k，L₀，l₀，L，N与所要模拟的高斯光束特征以及大气条件和传输距离有关，根据不同的需要选取不同的值.为便于比较，本文采用与文献[20]中相同的参数设置，即考虑某型地基外差激光雷达在近地面工作时分别遭遇强、弱湍流的情况^[20]，具体参数为λ=2 μm，$C_n^2$=10^-12 m^-2/3, Z=0~2 000 m($C_n^2$=10^-15m^-2/3, Z=0~10 000 m)，L₀=3 m，l₀=1 cm.根据采样约束条件^[7]，本文选择相位屏采样间隔为0.002 m，采样点数N=512，则相位屏尺寸为L=1.024 m.

首先根据上述参数，利用上一节的模型生成高精度的相位屏，并对其结构函数进行统计计算与理论结构函数进行比较，经优化得到的四个参数值如表 1所示.

文献[1]中提到，利用该方法生成的相位屏精度不随湍流强度的变化而改变，因此表 1中的参数可用于不同强度的湍流.将表中参数代入上节的相位屏模型中，即可生成高精度的相位屏，图 4为由该模型生成的相位屏的一个实现，可以看见代表低频成分的相位屏倾斜较为明显，说明低频成分得到了补偿.

图 4. 相位屏生成结果

Fig. 4. A realization of simulated phase screen

下载图片 查看所有图片

为考察所生成的相位屏精度，考虑弱湍流条件下$C_n^2$=10^-15 m^-2/3，2 μm波长的光束传输2 000 m的情形，利用相位屏模型独立生成了1 000幅相位屏，统计计算其结构函数并与理论结果进行比较，如图 5所示，作为参照，利用FFT方法、次谐波补偿方法以及改进前的最优化方法也生成了1 000幅相位屏，统计计算其结构函数.

图 5. 三种方法生成的相位屏的结构函数与理论结构函数对比

Fig. 5. The structure functions of phase screens generating by different method vs. theoretical structure functions

下载图片 查看所有图片

从图 5可以看出，利用改进最优化方法生成的相位屏，其结构函数与理论结果十分符合，作为对比，FFT方法生成的相位屏除了在开始的小尺度范围内与理论结果符合很好外，在其他尺度范围上与理论结构函数偏离很大；经过次谐波补偿后的相位屏在0~0.25 m范围内与理论结果符合很好，但是超过这一范围后，其结构函数将会超过理论结果，造成过度补偿，这是由于次谐波方法仅能校正零频附近的低频误差，校正能力有限，无法预测结构函数在大尺度上的变化速率衰减；改进前的最优化方法具有一定的误差校正能力，但校正程度不如改进后的方法且仍存在过度补偿的问题.经计算，改进最优化方法生成的相位屏，其结构函数与理论结果在大尺度范围内的最大相对误差为0.93%，而改进前的最优化方法得到的结果为6.75%(优化参数c₁=7.430, c₂=51.628)，FFT方法为22.99%，次谐波法为16.81%.可见我们提出的最优化方法与传统方法相比能够有效降低相位屏的结构函数误差，并且利用改进后的最优化方法还能够有效纠正其他方法产生的过度补偿问题.

利用所生成的相位屏，对上文提到的某型地基外差激光雷达在不同大气条件下的传输进行了模拟，其中激光波长为2 μm，腰斑半径为W₀=0.07 m.模拟算法采用split-step方法^[6]，即将传输路径上的大气分割成不同的湍流层，每个湍流层以相位屏表示，激光在该层的传输分为真空传播以及相位屏调制两个独立过程，每层湍流所代表的Rytov方差小于整个路径Rytov方差的0.1倍^[21]，按照这个要求，针对$C_n^2$=10^-15 m^-2/3的情况，我们将0~10 000 m路径上的湍流分为20层；针对$C_n^2$=10^-12 m^-2/3的情况，将0~1 000 m路径上的湍流分为20层，将1 000~2 000 m路径上的湍流分为32层.图 6所示为波长为2 μm的激光束分别在$C_n^2$=10^-15 m^-2/3和$C_n^2$=10^-12 m^-2/3时传输1 000 m后的光强分布，分别代表弱湍流情形和强湍流情形.

图 6. 弱湍流和强湍流情况下的光斑强度分布

Fig. 6. Beam irradiance distributions in weak and strong fluctuation levels

下载图片 查看所有图片

从图 6中可以看出，激光束在弱湍流情况下光斑形状较初始的高斯形略有不同(白色椭圆表示未经大气湍流扰动的激光束轮廓)，光斑瞬时中心发生了偏移(白色圆点为畸变后光斑的瞬时光强中心位置)，但是并未发生光斑破碎；而在强湍流情况下，光斑碎裂成很多子光斑，此时光斑的瞬时中心已无明显的偏向性.这与理论分析是一致的，光斑瞬时中心的偏移在弱湍流情况下较为明显，而随着湍流强度的增加，光斑偏移情况越来越弱，逐渐趋于饱和.但是单纯的光斑分布无法定量的描述光束经大气湍流后的统计特性，因此需根据模拟结果计算激光束的光强统计量.本文选择两种二阶统计量进行分析，包括光束扩展(Beam spread)、光束漂移(Beam wander), 为减少光强的不确定性，对400个样本进行统计，可将不确定性控制在1%以内^[22].

3 光束扩展及光束漂移计算结果与分析

激光束在湍流大气中的光束扩展使得激光中心能量发生弥散，导致中心能量降低，将会影响激光通信、激光雷达等系统的工作效率^[23-24].光束扩展受光束在真空中的衍射扩展、小尺度湍涡对光束的额外衍射扩展以及大尺度湍涡对光束造成的折射三种因素的影响，这里的光束扩展是在长曝光图像中观测到的，所以也叫长期光束扩展^[25]，而光束在真空中的衍射扩展以及小尺度湍涡引起的额外扩展被称为短期光束扩展^[26]，大尺度湍涡的折射效应引起的光束扩展被称为光束漂移^[27].长期光束扩展与上述三种因素之间的关系可表示为^[28]

式中，$W = {W_0}\sqrt {\Theta _0^2 + \Lambda _0^2} $为高斯光束在真空中传输Z距离后的光斑半径，${\Theta _0} = 1 - \frac{Z}{{{F_0}}}$和${\Lambda _0} = \frac{{2Z}}{{kW_0^2}}$被称为输入平面光束参数，F₀为高斯光束在输入平面处的曲率半径，W₀为腰斑半径，T_SS表示小尺度湍涡的贡献，T_LS表示大尺度湍涡的贡献.在解析表达式的推导上，ANDREWS L C等^[25]利用广义惠更斯-菲涅尔原理得到了适用于弱湍流条件下的长期光束半径，表示为

式中，$\sigma _{\rm R}^2 = 1.23C_{\rm n}^2{k^{7/6}}{Z^{11/6}}$为Rytov方差，$\Lambda = \frac{{{\Lambda _0}}}{{\Theta _0^2 + \Lambda _0^2}}$与另一参数$\Theta = \frac{{{\Theta _0}}}{{\Lambda _0^2 + \Theta _0^2}}$被称为输出平面光束参数，从式(5)可以看出$T = {T_{\rm SS}} + {T_{\rm LS}} = 1.63\sigma _{\rm R}^{12/5}\Lambda $.光束漂移方差可利用有效光束参数法^[29]推导得出，其结果为

式中，${\overline \Theta _0} = 1 - {\Theta _0}$，式(6)可用数值积分的方法以很高的精度进行求解，则光束漂移引起的扩展为$\sqrt {2\left\langle {r_c^2} \right\rangle } $.短期光束扩展可表示为${W_{\rm ST}} = W\sqrt {1 + {T_{\rm SS}}} $，根据式(4)所示的关系可直接得到W_ST的表达式为

至此我们已得到与光束扩展相关的所有量的解析式，而对于从模拟结果统计计算上述三个参数的问题，可通过求解光强的二阶矩得到^[30]

式中，

为光强的重心位置即瞬时重心位置，|r|²=x²+y²，假设〈r_c〉=0，则$\left\langle {r_{\rm c}^2} \right\rangle = \left\langle {r_{{\rm c}x}^2} \right\rangle + \left\langle {r_{{\rm c}y}^2} \right\rangle $即为光束漂移方差.而长期光束扩展则可表示为$\sqrt {2{{\left\langle {{\sigma ^2}} \right\rangle }_{{r_{\rm c}} = 0}}} $，短期光束扩展则可表示为$\sqrt {2\left\langle {{\sigma ^2}} \right\rangle } $.

利用式(8)和式(9)，对模拟得到的两种大气条件下的光束进行了统计计算，并与式(4)~(7)所示的理论结果进行比较，同时为比较修正大气谱与Von-Karman谱的不同，还引入了Von-Karman谱的模拟结果，如图 7至图 9所示.

图 7. 长期光束扩展模拟结果

Fig. 7. Simulated long-term beam width

下载图片 查看所有图片

图 8. 短期光束扩展模拟结果

Fig. 8. Simulated short-term beam width

下载图片 查看所有图片

图 9. 光束漂移模拟结果

Fig. 9. Simulated beam-wander standard deviation

下载图片 查看所有图片

由图 7可知，对于长期光束扩展来说，在弱湍流至中等湍流情况下(图 7(a)所示的弱扰动情况以及图 7(b)所示的强扰动情况下传输距离较短时)，模拟结果与理论结果符合得很好，但是当湍流强度由中等湍流至强湍流变化时(图 7(b)所示的强扰动情况下传输距离较长时)，模拟结果将小于理论结果，并且随着传输距离的增加，模拟结果与理论结果偏差越来越大.其原因在于式(5)所示的长期光束扩展解析式的推导是以高斯光束经过湍流过后仍近似为高斯光束为前提的，这一假设在弱湍流情况甚至中等强度湍流情况下都是成立的^[31]，但是当进入强湍流区域时，由于光斑碎裂为许多的小光斑，使得原来的高斯型光束变成非高斯型，于是这一假设便失效了.而通过对图中采用修正大气谱与采用Von-Karman谱所得的模拟结果比较可知，在弱湍流情况下，采用两种谱所得的模拟结果相差不大，而在中等湍流以及强湍流情况下，采用修正大气谱所得的模拟结果大于采用Von-Karman谱的结果，并且在中等湍流中与理论结果更相符.

图 8所示的短期光束扩展模拟结果与图 7所示的结果非常相似，因此也与图 7的分析过程是一致的.而图 9所示的光束漂移的模拟结果与光束扩展的情况有所不同，由图 9可看出，在弱湍流(图 9(a))以及中等湍流(图 9(b)中传输距离小于1 000 m)情况下，采用修正大气谱所得的模拟结果大于采用Von Karman谱所得的模拟结果，并且与理论结果符合得更好；而当图 9(b)中传输距离大于1 000 m时，光束漂移方差的增长速度将随着传输距离的增加变得越来越慢，并最终趋于饱和，这是由于光斑的破碎使得光束漂移对湍流效应的贡献越来越小.从图 9还可看到，当采用修正大气谱时，光束漂移方差趋于饱和的速度快于Von Karman谱，这是因为修正大气谱的“凸起”加快了光斑的破碎.

从上述分析中可得，以上三种统计量在弱湍流至中等湍流情况下与理论结果符合得较好，而在强湍流情况下，模拟得到的结果与理论结果发生偏离，并且小于理论预测的结果，这种偏离可用光斑在强湍流中的破碎现象进行解释.在与Von-Karman谱的模拟结果进行比较时发现，本文所得的模拟结果略大于采用Von-Karman谱的结果，这与修正大气谱存在的“凸起”直接相关，也印证了修正大气谱的“凸起”会反映在光束扩展的预测上，而利用本文提出方法生成的相位屏能够很好地预测这一特征.

4 结论

修正大气谱在高波数处有“凸起”的存在，能够很好地表示实际大气折射率扰动特点，因此以该谱为基础生成的相位屏也将更好地反映出光束经湍流效应后所具有的各种特征.本文提出了一种高精度的相位屏生成方法，该方法以文献[1]中的基于优化的相位屏模型为基础，通过对低频区域的扩展以及采样点重设置，将原有的单一补偿屏扩展为双补偿屏，然后利用引力搜索算法寻找低频采样区的最优采样点分布，进而克服低频采样不足的缺点.仿真结果表明，本文提出的相位屏生成方法可将低频区最大相对误差降低至1%以下，而经典FFT方法低频误差可达23%，次谐波法为16.81%，即使是相同的算法，在改进采样设置之前也有6.75%的相对误差.同时，利用本文方法生成的相位屏对某地基外差激光雷达在两种湍流条件下的传输进行了模拟，并计算了其二阶统计量(光束扩展以及光束漂移)，结果表明，在弱湍流至中等湍流范围内，模拟结果与理论结果十分一致，而在强湍流情况下由于光斑破碎使得模拟结果会低于理论结果的预测，这是由于理论公式是在未考虑光斑破碎的情况下推导出来的，在强湍流情况下将会失效.在与Von-Karman谱的模拟结果比较时可以清晰地看到由于修正大气谱存在“凸起”，使得最终的统计结果略大于Von-Karman谱的统计结果，说明本文提出的相位屏生成方法能够有效的表征实际大气折射率的扰动.