Simulation method of pulse laser fuze echo in dust environment

大气中的烟尘对以激光为探测手段的各类武器系统的性能有着不可忽略的影响。特别是在复杂的战场环境下，由于爆炸或燃烧引起的烟雾、沙尘粒子会对激光信号产生折射、吸收、反射等作用，从而使得激光束能量发生衰减和空间分布的改变。尤其在激光近距探测中，由于大气中悬浮烟尘粒子的后向散射产生的回波干扰信号，使得激光探测设备产生误判和虚警，影响了引信系统的性能和精度，这一问题已成为限制激光近距探测发展和应用的瓶颈问题^[1-5]。且微粒的散射特性随着自身物理特性及激光波长改变，表现出特有的性质，因此，微粒的光散射研究对于各种激光系统具有重要的实际价值，尤其在设计激光近距探测系统时，必须要考虑激光在大气中的传输及衰减。

研究大气中烟尘环境传输特性的主要方法有渐进法、离散纵坐标方法、蒙特卡洛方法等，蒙特卡洛方法基于辐射传输方法但无需求解辐射传输方程，而是直接模拟光子的传输过程，可高精度地处理多种类型的辐射传输问题，因此常用在微粒的多重散射研究中，而蒙特卡洛方法离不开粒子散射相函数的模拟计算。目前在蒙特卡洛仿真时通常利用H-G散射相函数进行散射角抽样，而H-G散射相函数仅与粒子不对称因子有关，只能近似描述粒子的散射角分布，尤其是后向散射部分与实际差异较大，而T矩阵法可对粒子散射特性进行精确计算，但在蒙特卡洛仿真研究时，大部分学者仍采用数值计算简易的H-G函数，也未给出H-G散射相函数的精度及适用范围；且在激光烟尘环境多重散射研究中，各种理论分析占主导位置，相应的实验研究及验证较少。因此文中基于T矩阵法及半解析蒙特卡洛仿真，提出一种基于离散散射相函数的散射角抽样法，建立激光引信在烟尘环境发射及接收的理论模型，仿真分析不同烟尘浓度及湿度环境对激光回波的影响，并进行不同浓度烟尘环境的激光引信回波实验。一方面通过实验分析实际工作中不同烟尘浓度对激光引信性能的影响，另一方面将数值模拟结果与实验结果进行对比，验证了文中采用的T矩阵散射相函数分析激光多重散射的准确性，为复杂烟尘环境的激光引信实验技术及精确蒙特卡洛数值模拟奠定理论基础。

常用的散射相函数计算方法有：Legendre多项式展开法、H-G散射相函数、Mie散射理论、T矩阵等，而在蒙特卡洛方法中，为方便抽样光子碰撞后得到的散射角，大都采用表达形式简单、数值计算方便的H-G散射相函数，表达式如下：

式中：$\theta $为散射角；$g$为不对称因子。

T矩阵法不受入射场和散射场的约束，只与散射粒子的形状、大小、折射率及坐标系中的位置有关，即T矩阵只需计算一次，就可以得到任意入射波产生的散射场^[6-7]。根据麦克斯韦方程的线性以及本质关系，散射场展开系数${p_{mn}}$、${q_{mn}}$和入射场展开系数${a_{mn}}$、${b_{mn}}$之间的关系也必然是线性的，可由转换矩阵（T矩阵）给出：

假设所求粒子在某个取向时的T矩阵已知，根据Wigner D函数的正交性，且子矩阵${{{T}}^{ij}}$(i, j =1,2)都是对角阵，由此可以方便地计算出衰减截面为：

和散射截面为：

则粒子群的衰减系数及散射系数可写为：

式中：$N$为粒子的数量浓度。且根据广义球函数的正交性，可得到粒子的不对称因子为：

式中：${a_1}(\theta )$即为利用勒让德多项式展开的散射相函数，表示为：

式中：$g_{00}^s$、$g_{0 - 0}^s$均为扩展系数，对于旋转对称粒子，散射相函数满足以下归一化条件：

利用上述两种方法分别计算几种粒径烟尘粒子的散射相函数，如图1所示，其中相应的粒子散射强度分布的极坐标图也展示出，可以更直观地观察不同粒径粒子光散射的角分布。由于H-G函数的局限性，计算中假设粒子形状为球形。烟尘粒子半径范围为0.1~10 μm^[8]，激光波长选取较为常见的905 nm，干燥粒子复折射率$m = 1.520 + 0.008i$^[9]。由图1可知，随着粒子粒径的增加，T矩阵法计算得到的散射相函数随着角度的变化波动逐渐剧烈，而H-G散射相函数比较平滑。且随着粒子粒径的增加，两者在小角度（0°~90°）及大角度（150°~180°）附近散射强度的差距也逐渐增大，因此在蒙特卡洛仿真粒子散射角抽样时，H-G散射相函数无论是在透过率亦或是后向散射模拟上，均会与T矩阵散射相函数的抽样角度有很大差异。

Figure 1.  Comparison of two scattering phase functions and scattering intensities. (a) Particle radius: 0.1 μm; (b) Particle radius: 0.4 μm; (c) Particle radius: 0.8 μm; (d) Particle radius: 2 μm; (e) particle radius: 6 μm; (f) Particle radius: 10 μm

蒙特卡洛数值模拟本质上是通过相应的概率模型和随机数来模拟单个光子在介质中的随机行走过程，并通过对大量光子的追踪得到光在介质中传输行为的统计结果^[10]。由于T矩阵得到的是离散相函数，无法像H-G函数那样求解出散射角的逆函数，因此提出基于T矩阵散射相函数的散射角抽样方法，使其能很好地运用到蒙特卡洛方法中。

在蒙特卡洛数值计算中，散射角必须满足概率密度分布函数，对于常用的连续散射相函数H-G函数，可通过公式(10)由0~1之间的随机变量$\xi $得到散射角为^[11]：

T矩阵法较难获得对应散射相函数的解析解，无法像H-G函数那样解析积分求逆，因此对于T矩阵得到的散射相函数，在散射角间隔范围内取n个分割点（文中共取1 801个分割点），即可得到第i个分割点处的散射相函数值${a_1}({\theta _i})$，令：

借助概率基本定律，同样利用均匀的随机函数${\xi _1} \in (0,1)$，建立${\xi _1}$与$i$两个随机变量之间的对应关系，即寻找未知数$m$得到$\min \left| {{P_m}(\theta ) - {\xi _1}} \right|$，即可得到仿真所需的随机散射角${\theta _m}$及对应的散射相函数${a_1}({\theta _m})$。假设粒子此时的散射方位为$(u,v,w)$，碰撞后，通过抽样得到散射角${\theta _m}$，粒子新的散射方位可由下式表示为：

但当散射角很小时，即$\left| w \right| > 0.999\;99$，新的散射方位就需使用下式：

图2对比了图1(c)中T矩阵法计算的散射相函数抽样角度分布与H-G散射相函数抽样角度对比，与图1(c)结果一致，在小角度及（0°~90°）及大角度（150°~180°）附近差异较大，且从结果中可分析出，采用H-G相函数模拟透过率会得到较大的结果，而在回波分析时仿真数值会小于T矩阵散射相函数得到的回波。

Figure 2.  Comparison of two sample methods of scattering phase functions

将激光在烟尘中传输的问题转化为大量光子在烟尘粒子中的输运问题，光子的每一次散射只与前一次散射状态有关，统计光子经过烟尘粒子散射和吸收后，激光接收系统探测到的光能量。

建立如图3所示笛卡尔坐标系：发射和接收视场轴线平行，为方便计算，以探测器接收中心为坐标原点，激光接收与发射连线为x轴，光束发射方向为z轴。定义发射中心与接收中心距离为${d_0}$，发射光束半径为${r_e}$，发散半角为${\varphi _e}$，激光接收口径为${r_r}$，视场半角为${\varphi _r}$。由于激光发射口径较小，因此将发射光源假设为点光源，每个光子初始位置即为坐标$ （{d}_{0},0,0）$，初始散射角${\theta _0}$在$\left[ {0,{\varphi _e}} \right]$内均匀分布，初始方位角${\phi _0}$则在$\left[ {0,2\pi } \right]$内均匀分布。将发射脉冲在时间上以$\Delta t$为间隔划分为光子包，发射光脉冲峰值光子数为${N_0}$，则每个光子包中光子的数量表示该时刻的发射功率。

Figure 3.  Emission and receive model in dust environment

光子在烟尘环境中传输时，可能发生吸收或散射两种过程，利用单次反照率${\omega _0}$与随机数${\xi _2} \in (0,1)$进行判断，其中${\omega _0} = {\mu _s}/{\mu _t}$。若${\xi _2} > {\omega _0}$，光子被吸收，反之则发生散射。光子在两次碰撞间的随机步长遵循朗伯-比尔定律，可利用随机数$ {\xi }_{3}\in （0,1）$表达，即$\Delta s = {{ - \ln {\xi _3}} / {{\mu _t}}}$。在多重散射模拟中，光子的每一次散射只与前一次散射有关，发生散射后的运动方向即由散射相函数抽样确定，而每次的散射方位角都可认为在$(0,2\pi )$均匀分布。

通过上述步骤，可以对光子路径进行追踪，得到光子每次散射后的坐标位置，而在回波仿真中，最关心的是能被探测器接收的光子。为提高计算效率，采用半解析接收法^[12]，只要光子位于接收视场且权重高于阈值${\omega _e}$，则进行一次半解析接收，建立如图4所示的半解析接收几何模型。

Figure 4.  Semianalytic sensing geometric model

则光子每次返回接收系统的概率可以表示为：

式中：$\theta '$为接收时光子散射方向与接收方向的夹角；${a_1}(\theta ')$为该夹角下对应的散射相函数；$\Delta \Omega = A/R_m^2$，为接收时的立体角，其中$A$为激光接收系统面积，${R_m}$为光子散射位置距接收系统中心的直线距离；${{\rm e}^{( - {\mu _t}{R_m})}}$为光子直接运动${R_m}$而不发生碰撞的概率；${\omega _m}$为接收时光子的权重，光子进行一次半解析接收后，新的权重不再采用${\omega _n} = {\omega _{n - 1}}{\omega _0}$，而是变为：

半解析接收后，只要光子权重仍大于阈值权重，则继续追踪光子。

每次在新的散射发生前，除了对光子的吸收及散射过程进行判决，还需进行光子终止条件的判决。(1)若光子离开烟尘边界，则停止追踪；(2)光子权重低于阈值${\omega _e}$，则判定光子消亡；由于采用半解析接收方法，因此需要增加光子的终止条件，即(3)若光子连续两次散射处于接收视场外，则判定光子逃逸，开始追踪下一个光子。持续追踪每个光子的运动轨迹，根据接收时光子在烟尘中总的传输路径进行分段，统计每段的概率之和，最终通过对大量光子的模拟得到接收器接收到的信号波形。

为分析不同浓度烟尘对激光引信传输性能的影响，并验证理论模型计算激光回波的准确性，建立烟尘环境激光引信实验室，如图5所示。

Figure 5.  Dust environment experiment. (a) Laboratory layout; (b) Experimental scene

实验室主要分为三大部分：控制室、测试区及沉降室。控制室用于设备安置及调试，方便实验人员记录数据；测试区产生实验所需的烟尘环境，且测试区距离可调节；沉降室便于烟尘的清理。实验室长度为20 m，宽度及高度均为2 m，可视性和密闭性较好。利用发烟饼（主要成分为铵盐），制造无风烟尘环境，并记录烟尘实时浓度及激光接收回波。实验中采用的脉冲激光发射波长为905 nm，峰值功率75 W，实时粉尘监测仪为CASELLA公司的CEL-712 Microdust Pro，最大量程为250 g/m³。

为了更好地对比实验与理论模型计算结果，仿真参量与实验系统设置一致，具体如表1所示，其中，$\tau $为发射激光脉宽。

利用CEL-712 Microdust Pro配套的Insight数据管理软件处理烟尘数据，由于烟尘浓度实时变化，因此取相对稳定时段的平均值（去除极值）作为某一时刻的烟尘浓度，并对应同一时间的激光回波信号作为一次实验结果。此次实验主要研究烟尘浓度变化对激光回波的影响，因此在测试区尽头使用吸光材料，模拟无目标情况下烟尘环境造成的回波，且实验时引信处于烟尘环境中，因此实验结果与距离无关。实验结果如表2所示，随着烟尘浓度的增加，激光回波幅值也相应增大。

假设粒子成分单一，根据封闭环境铵盐烟尘实验经验公式^[8]可得质量浓度$C$与数量浓度$N$间的关系为：

以及

式中：${k_0}$为凝结常数；$\alpha $、$\;\beta $均为损耗常数。文中取对应无风环境时$\alpha = 3.4 \times {10^{ - 3}}$，$\beta = 3.2 \times {10^{ - 3}}$，${k_0} = 2.2 \times $$ {10^{ - 8}}\;{\rm{c}}{{\rm{m}}^{ - 3}}$，初始值取测量值${C_0} = 1.02\;{\rm{g/}}{{\rm{m}}^{\rm{3}}}$，${N_0} = 2.5 \times $$ {10^7}\;{\rm{c}}{{\rm{m}}^{ - 3}}$。则可得相应质量浓度下的数量浓度，如表3所示。

为进一步分析湿度对烟尘环境激光引信回波的影响，采用Hänel总结的粒子粒径及复折射率随相对湿度变化的规律，如下式所示^[13]：

式中：${m_r}$、${m_i}$为干燥粒子复折射率的实部和虚部；${m_{rw}}$、${m_{iw}}$为水复折射率的实部和虚部；${m_{re}}$、${m_{ie}}$为湿粒子复折射率的实部和虚部；$r$、${r_e}$分别为干燥粒子半径和湿粒子等效半径；$frh$为吸湿增长因子；${\rm{RH}}$为相对湿度；$\upsilon $为质量增长粒子，文中取值3.9。文中选取的入射波长为905 nm，此时纯水的复折射率为1.328+5.12×10⁻⁷i，因此粒子在高湿度（RH=95%）情况下复折射率为1.439+4.37×10⁻³i、中湿度（RH=50%）情况下复折射率为1.3458+6.65×10⁻⁴i。

分别计算高湿度（RH=95%）、中湿度（RH=50%）以及干燥情况（RH=0%）下烟尘粒子的散射相函数，如图6所示，可知，当环境湿度发生变化，由于粒子的吸湿性，其散射特性也发生变化，虽然整体趋势比较接近，但随着相对湿度增加，前向散射及后向散射均增强，尤其是高湿度情况下更为明显，因此，在采用回波信号的脉冲激光近距探测研究中，环境湿度在一定程度上会影响探测结果。

Figure 6.  Scattering phase functions with different relative humidity conditions

根据计算得到的散射相函数利用文中所述散射模型仿真不同烟尘浓度在三种湿度环境下的激光回波信号，如图7所示，图为仿真得到的反映回波功率变化的回波波形，为方便比较，按回波最大值归一化处理。

Figure 7.  Normalized simulation echo powers with different relative humidity conditions

结果显示，随着湿度增加，烟尘回波也明显增大，这是由于湿度增加时，粒子吸湿因子增大，其散射特性也相应增强，使得探测器接收到的光子数增加。

而随着烟尘浓度增加，三种湿度条件下，激光回波均增大，这是因为粒子密度增加，后向散射的次数也相应增多，能到达接收视场的光子总数就增加，与实验结果一致，这一点并不因为湿度变化就改变，如图8所示。且随着烟尘浓度的增加，湿度变化造成的回波增量也不断增大，因此，在烟尘浓度较小（<150 mg/m³）且对精度要求不是很高（<0.2）时，可不考虑湿度对激光探测的影响，但随着烟尘浓度的增加，环境湿度的影响因子也会增加，因此在仿真及实验时都需关注。

Figure 8.  Normalized simulation echo powers changing with dust concentrations in different relative humidity conditions

实验时间为冬季，南京室内约为50%中湿度条件，为验证模型的准确性，利用前文所述模型并分别采用两种散射角抽样方法仿真对应烟尘浓度下激光回波信号，结果如图9所示，同样为方便比较，按T矩阵法最大回波进行归一化处理。

Figure 9.  Normalized simulation echo powers changing with dust concentrations. (a) T-matrix; (b) H-G scattering phase function

由图9可知，两种方法得到的数值趋势与实验结果一致，均随着烟尘浓度的增加而逐渐增大，且随着浓度的增加，波形脉宽也加大，这是由于随着光子与烟尘粒子碰撞次数增加，到达接收视场的时间延迟也越明显。而通过两种方法仿真结果对比可知，采用H-G散射相函数抽样的方法结果远小于T矩阵散射相函数抽样的结果，通过上文散射相函数对比可知，H-G散射相函数大角度（150°~180°）附近散射数值较小，而在回波仿真中，由于是多重散射，不仅是后向散射（180°）光子对回波有贡献，大角度散射处光子均有较大概率能到达探测器，仿真结果与理论分析一致。

由于实验得到的波形是经过光电转化的电压信号，且受光学系统透过率、探测器偏置电压等复杂因素影响，因此，仿真波形难以与其直接对比，但两者趋势应一致，且成比例关系。为进一步验证理论模型的准确性，将文中仿真结果与实验结果及常用的H-G散射相函数仿真得到的数值按各自均值归一化处理并进行对比，如图10所示，并计算仿真与实验相对误差，结果如表4所示。

Figure 10.  Comparison of results of two simulation methods and experiment

可知当浓度较小（<138.6 mg/m³）时，H-G散射相函数可得到较准确的仿真结果，而T矩阵散射相函数得到的结果相对误差较大，这是由于烟尘浓度较小时，实验烟尘环境会较不稳定，粒子流动性较大造成浓度不均匀，因此粒子后向散射次数反而减少。而随着烟尘浓度增加，实验室烟尘趋于稳定，T矩阵散射相函数的相对误差减小（<10%），而H-G散射相函数仿真得到的相对误差基本随着浓度增加而增大，这是由于H-G散射相函数仅是散射相函数的近似模拟，难以真实反映粒子的散射特性，特别是大角度附近散射误差较大，因此在计算回波时有一定劣势，浓度越大越差异越明显，而T矩阵散射相函数更为精确，得到的回波幅值与实验结果较接近，因此总体来说基于T矩阵方法得到的散射相函数在激光回波研究中有更高的精确度。

基于T矩阵法计算不同粒径烟尘粒子的散射特性，包括散射系数、消光系数、不对称因子、散射相函数等，结合烟尘环境下激光引信发射及接收模型及散射相函数的物理特性，提出了基于离散散射相函数的光子散射角度抽样方法。建立激光引信在烟尘环境中的回波仿真模型，利用激光引信近距探测和接收视场小的特点，简化了散射光子在空间的概率分布，仿真了不同浓度不同湿度条件的烟尘环境下905 nm脉冲激光回波信号，并对仿真结果进行了实验验证。结果表明，烟尘粒子的后向散射特性随湿度的增加而加强，继而使得回波信号增大，通过对比几种浓度及相对湿度下的回波信号可知，当烟尘浓度较大（≥150 mg/m³）、精度要求较高（≤20%）时，必须考虑湿度变化对激光回波的影响；同时，在烟尘环境激光回波仿真中，传统的H-G散射相函数操作简便，且在烟尘浓度较小（<138.6 mg/m³）时精度较高，但随着烟尘浓度增加误差也逐渐增大，因此在烟尘环境单一且烟尘浓度较小的情况下，仍可采用H-G散射相函数进行多重散射的蒙特卡洛仿真；而采用T矩阵散射相函数仿真的结果与实验吻合度较好，尤其在大浓度烟尘环境下，相对误差均小于10%，与H-G散射相函数相比精度较高，因此基于T矩阵散射相函数的蒙特卡洛方法模拟烟尘环境下的激光回波信号是可行及准确的。文中研究可进一步用于复杂烟尘环境下激光引信的性能分析及优化设计。