中高层电参数的昼夜差异对地闪回击电磁场传播的影响

杨许铂，黄翰林

（1.海南热带海洋学院 海洋信息工程学院，海南 三亚 570022；2.沈阳理工大学 理学院，沈阳 110159）

0 引言

因为大电流和微秒量级的短暂时间尺度，云地闪回击辐射出强烈的电磁脉冲（electromagnetic pulse，EMP）。该辐射场能量主要集中在甚低频（Very low frequency，VLF）波段[1］，可长距离传播。一部分辐射波向上传播，与高层大气中电离层区域的等离子体相互作用，加热此处的电子，增加电子与中性分子的碰撞率，增强其电离和吸附作用，并激发光化学反应，产生红环精灵（Elves）。同时，闪电会中和雷暴云中的电荷，打破静电平衡，产生准静电场（Quasi-electrostatic field，QE），QE 感应场同样可渗透进入电离层底层，加热引起电离和光化学反应，产生红色精灵（Sprits）[2］。QE 和EMP 的加热机制也许相同，但二者的场形态特征和持续时间尺度不同，因此对电离层的加热和电离区域有所不同，产生的Elves和Sprits也不一样。

闪电放电电场作用于电离层低层是一个复杂且具有多特征的物理问题，对其进行科学探究时除了观测实验外，还需要借助数值模拟方法以探讨其内在的物理机制图像。Taranenko 等[3-4］首次建立了一个一维全动力学模型模拟EMP 场引起的电离作用，计算了其光辐射含量变化。Inan 等[5］则建立一个二维圆柱坐标模型，采用时域有限差分方法来模拟分析地闪EMP 引发的电离层参数变化。Veronis 等[6］、Barringtonleigh[7］扩展了这个二维模式，考虑约500 ms 特征时间的慢速放电电流，研究了感应场对电离层的影响。Cho等[8］分别编写了一个电磁场和一个准静电场的计算机程序，考察辐射场和感应场的结构形态以及产生的Elves和Sprites光学特征。Rodger等[9］采用位于美国检测到的闪电数据，带入Cho等[8］的电磁场程序中，模拟研究发现连续的闪电EMP 事件可能会使得低电离层电子密度增大至背景值2 倍以上。Marshall 等[10］提出了一个三维的闪电EMP 与低电离层相互作用的时域有限差分模型，该模型预测连续的云地闪和云闪活动会引起电离层D层密度的长时变化。张其林等[11］则计入了强辐射场中高层大气的非线性加热与电离作用，使用时域有限差分模型模拟辐射场在中高层大气的波形，并和Lu[12］的解析算法进行了对比分析。Wang等[13］采用张其林等[11］的模型，计入重力波对中性波大气的扰动，模拟研究该条件下辐射场波形的变化及光学特征。

上述模型较好地解释了一些闪电放电作用于电离层效应的观测现象，但这些模型使用的是简单恒定大气电导率剖面以及仅两个参数确定的双指数函数的电子密度公式[14］，没有考虑时间变化和地区差异，一定程度上限制了模型的研究范围。根据麦克斯韦理论，大气中电导率决定电磁波能量损失的程度，其剖面形态更决定着电波反射高度，电子密度也是非线性效应中的重要参数。国际参考电离层（IRI）经验模型可在指定时间地点参数下提供电离层电子离子及温度等信息，Schlegel等[15］给出了一种已知电子密度及温度等参数下计算D 层电导率的方法，本次研究中将使用该方法获得可随时间地点变换的大气电导率剖面及电子密度剖面。

从麦克斯韦方程组出发，采用时域有限差分法构建一个二维地闪回击激发的电磁波在地面-电离层波导传播模型，模型耦合了高层大气电场对电子加热和电离的非线性效应。其背景参数则采用基于国际参考电离模型提供的电离层电子等参数以及用这些参数构建的大气电导率剖面，并以此分析昼、夜条件下，地闪回击电磁场在电离层传播的差异。

1 地闪回击电磁波传播模型

1.1 模型出发方程

随着海拔高度增加，中性大气的带电粒子成分逐渐增多，地闪回击激发的电磁波在此区域传播时，需要考虑到导电大气对电磁波的损耗。可列下面麦克斯韦方程组[16］为

其中：σ为大气电导率；分别为电场场强和磁感应强度；ε0和μ0分别是真空中的介电常数和磁导率;为回击电流密度，竖直向上为正，位于海拔高度0～10 km 的闪电通道中，是闪电放电EMP 的源。基底电流Ir(0，t)和回击电流密度Jr(r，t)之间有关系式

如图1所示，建立一个圆柱坐标系(r，φ，z)，令地闪回击通道位于z轴上，其他参数皆为轴对称。因电磁波的源是回击电流=Jr，产生的磁感应强度矢量方向为ϕ方向，据此展开Maxwell方程（1）、（2）式，得到

图1 建模示意图（a）坐标系（b）离散网格划分

其中：Er、Ez分别表示水平方向和竖直方向的电场分量；Bϕ则为磁感应强度的环向分量；c为真空中光速。

采用时域有限差分法离散方程（4）～（6）[11，17］，化简后可得到出发方程组

该离散方程对应着二维时域有限差分法中的横电场模式。其中：下标i、k分别表示柱坐标中水平距离r方向和竖直距离z方向的格点数；上标n表示迭代的时间步数；Δr、Δz分别为r、z方向的步长距离；Δt表示时间步长。(Er)ni+1/2，k即是Er(r=[i+ 1/2]Δr，z=[k]Δz，t=[n]Δt)的缩写表示，其他物理量以此类推。式中电场取值为网格边的中点，而磁场取值为网格正中间，磁场电场相隔半个距离步长。同样，磁场计算在半时间步上，而电场计算在整数时间步。

当i= 0时，r0= 0，式（9）右边的分式中因分式存在r0的因子而无法进行迭代计算，需要进行校对，在轴线上展开积分方程，考虑安培定律有

离散化得到

其中I即为式（3）定义的回击电流。

1.2 回击电场对电离层参数加热与电离的非线性效应

一方面，强烈的地闪回击产生的电磁场会加热冷等离子体、改变电离层参数，进而影响到局地的电离层电导率系数。另一方面，电离层电导率系数的改变又会影响闪电产生的电磁脉冲传播，两者互相作用。为了准确模拟雷暴激发电磁场在雷暴上空电离层处传播，模型里需要考虑自洽加热及电离过程[11］。

电子密度因碰撞电离和吸附而改变，方程式为

其中：ne表示电子密度；νi、νa分别为电子的碰撞电离系数和吸附系数，可看作是强电场的非线性函数。其中vi由下列式子[7］表示为

其中：N0表示海平面大气分子数密度；N表示大气分子数密度；系数a0=-624.68，a1= 249.60，a2= -32.878，a3= 1.454 6。

吸附系数νa[2］为

其中：系数为b0= -3 567.0；b1= 1 992.68；b2= -416.601；b3= 38.729；b4= -1.35113。

从大气中间层到电离层低层，忽略了磁场的作用，电子贡献的电导率σe可近似表示为标量

其中：e为电子电荷量；μe为电子迁移率。其折合电场的非线性函数为[2］

其中：x= log(E/N)；c0= 50.970；c1= 3.026；c2= 8.473 3 × 10-2。

模型将在每一个时间步上更新电离层高度处的电导率，并再代入迭代方程（8）、（9）、（11）中进行下一步场强计算。

1.3 背景电导率参数

在大气层低层（0～40 km），几乎不存在自由电子，电导率σ由正负离子所贡献，模型采用下列经验公式[18］

在电离层高度处，电导率原则上是受到地磁场影响的张量，分为彼德森电导率σp，霍尔电导率σH和场向电导率σ//，可由下列理论公式求得[15］。

其中：B表示地磁场强度；e为电子电荷量；下标e表示电子；下标i表示离子；ven、vin分别表示电子和离子对中性分子的碰撞频率；ωe、ωi分别表示电子和离子的磁回旋频率；me、mi分别表示电子质量和离子质量。在电离层D 层区域，因为σ//＞σH＞σp，而闪电回击引起的电磁波主要能量集中在VLF 频段，频谱能量在大于10 kHz 时则迅速衰减[1］。大部分能量的电磁波频率小于该处电子中性和离子中性碰撞频率，90 km 以下的电子回旋频率也小于电子碰撞频率，此时可忽略磁场的影响，电导率可表示成标量σ//。带电粒子的碰撞频率的计算公式[19］为

其中：Te表示电子温度；n[N2]表示氮气分子浓度；n[O2]表示氧气分子浓度；n[O]表示氧原子浓度。

计算出式（20）～（22）中所需的电离层电子浓度，在给定太阳和地磁活动参数、事件发生的时间与地点的情况下，电子温度和密度可以由国际参考电离层经验模式IRI-2020提供[20］，中性成分浓度则由美国海军实验室非相干散射雷达经验模式NRLMSISE-00 提供[21］，这两个经验模型可在NASA 社区协调建模中心（https：//ccmc.gsfc.nasa.gov/）在线计算。本文模拟选取雷暴发生地点为海南澄迈（地理经纬度N19.0°，E110.0°），日期为2018 年5 月9 日，以当地时12 时作为白天事例，00 时作为夜间事例。因IRI-2020 模型得到的电子密度在白天海拔高度下界为65 km，夜间为80 km，由式（21）计算得到的电导率存在下界。因此设由式（18）计算所得高度上界为40 km，中间大气部分采用对数线性插值方式来填补，其结果如图2所示，其中虚线部分是将式（21）中电导率里的电子贡献部分替换成无电场加热情况下的式（16）。图2 可看出，两种方法计算结果基本上一致。

图2 背景电导率参数

1.4 边界条件及电流源设置

模拟区域设置为水平半径350 km，竖直高度为110 km，水平距离步长Δr= 1 km，竖直步长dz= 1 km，时间步长dt= 10 ns。因为大地电导率远大于地表大气电导率，可假设地表为完美导体，即有下边界条件：

设置上边界和右边界为一阶Mur吸收边界条件[8，22］，上边界条件为

右边界条件为

左边界条件即为式（12）在轴线上展开，回击电流通道可设置为[0，10 km]。回击电流可以用传输线（TL）模型来描述，基底电流函数可选用Heidler等指数模型，其表达函数为

其中：τr表示电流的上升特征时间；τf表示电流的下降特征时间；I0为最大电流强度；位于t=τr时刻。

2 模拟结果与分析

设置强回击电流峰值I0= 220 kA，上升时间τr= 10 μs，下降时间τf= 100 μs，回击电流速度为v=1.7 × 105V/m强度，采用图2计算得到的背景电导率剖面，以此分析昼、夜不同情况下强地闪回击在电离层处引起的场强特征。本次模拟的初始的电离低层电子数密度与整个空间的中性分子数密度同样选取IRI-2020 和NRLMSIS-00 经验模式，因IRI 得到的电子密度存在下界，以电子密度指数模型[14］计算50 km 处的电子密度值，采用对数线性插值的方式[23］根据边界处的电子密度来填充50～65/80 km 空间区域上的变化。电子密度指数模型公式为

其中：β为电子浓度随高度的变化的梯度参数；h′为低电离层的参考高度。根据文献[23］里的实测数据，参数选取白天β= 0.3 km-1，h′= 70 km；夜间为β= 0.5 km-1，h′= 88 km。

模型运行时要在大气中高层更新其电子电离率、吸附率、密度、电导率，其计算空间日间设定为50 km以上，夜间设定为60 km以上，以下是模式的计算结果（见图3）。

图3 回击后不同时刻的空间电场分布

模拟很好地显现出地闪回击的强电磁脉冲传播特征。0.2 ms时，图3中可以分辨出强辐射场从准静电场中分离，随后向四周进行传播，在半径r方向上，辐射场随着传播距离增加以1/r比例逐渐减弱。而在海拔高度方向上，电离层高度上电导率较大，地闪回击辐射场则将在电离层高度处反射。在0.6 ms图像中可以看到首次强辐射场引起的反射天波波形，显弧光形状，在0.8 ms图像中亦可看见由天波传播至地表被地表反射的二次反射波形。水平半径50 km 内区域的电场为准静电场，与回击电流转移电荷有关，时间尺度与整个回击电流尺度有关，因此在模拟时间范围内变化不大。

图3中亦可看出强辐射场传播的昼夜区别，日间电离层反射高度比夜间电离层反射高度低，这是因为昼夜电导率差异引起的。因昼夜太阳照射不同，各层电离层的电子密度往往比晚上大1个数量级，其中最显著的是的D层（60～80 km），日落后此处电子密度很快耗尽，D层消失。因此白天在60 km处辐射场开始受到弱等离子体的作用，但到晚上，辐射场需要上到海拔高度80 km左右才开始受到电子对它的影响。

图4 为辐射场（r＞50 km）的最大场强随高度变化的剖面图，夜间辐射波场强Emax（80 km）/Emax（60 km）值在0.4，0.6，0.8，1.0 ms 时刻分别为0.91，0.96，0.91，0.84，而Emax（90 km）/Emax（60 km）的分别是0.28，0.18，0.12，0.08；由此可判定，夜间回击引起的强辐射场在海拔高度80～90 km之间和电离层冷等离子体相互作用，电磁波在此处发生折射反射。同样，白天辐射场强Emax（80 km）/Emax（60 km）为0.60，0.48，0.38，0.31；而Emax（85 km）/Emax（60 km）值为0.28，0.19，0.13，0.10；表明在白天回击引起的强辐射场与电离层冷等离子体作用在海拔高度85 km以下，其中大部分能量在80 km以下，电磁波在此处发生折射反射。

图4 地闪回击后最大辐射电场高度剖面图

电场加热带电大气引起的电子的电离与吸附系数变化可看成折合电场E/N的函数，而当vi=va时，此时电场大小记作Ek（图5）。在高于此阈值的稳定电场中，dne/dt＞0，因电离率与电子密度ne成正比，ne将呈指数式增长，形成电子雪崩过程，这一过程又称为空气击穿。Ek大小和气体浓度N有关，可采用经验公式Ek= 3.6 × 106N/N0来表示[2］，空气分子数密度越高则要求击穿电场值越大。一般而言，大气中性分子数密度是随着高度指数性下降，海拔高度越高Ek值越小，而强地闪回击电场可能在电离层高度上大于这个阈值，在短时间引起电离层电子密度的扰动，并产生各类的瞬态发光现象。

图5 吸附系数和电离系数变化

大气中性分子数密度日间与夜间差异较小，因此各高度击穿电场阈值日夜差异不大。但因日间电导率远高于夜间电导率，地闪回击的强电磁场在夜间比白天更容易渗透进入海拔高的区域（图4），因而夜间相较于白天更容易产生电子密度增强和各类瞬态发光现象。图6（a）为夜间Ek/E时空上的分布图，图上显示，大约在回击发生后0.3 ms，辐射场波峰传播至海拔高度80 km 处，此处空气较为稀薄，辐射电场强度足以发生击穿。但因击穿电场刚刚发生，电离和吸附引起的电子密度扰动需要有时长积累，因而此时电子密度几乎不变。0.4 ms，0.5 ms 时击穿电场区域分为2 部分，一部分随着辐射场往外扩散，另一部分位于雷击通道上空区域80 km 处，为准静电场渗透引起的击穿场。0.6 ms 至0.7 ms 准静电场因局地弛豫时间(τt=σ/ε)短暂而减小至击穿场强阈值下，辐射场则随传播距离r增大减少至击穿场强阈值附近，其击穿电场分布范围约为半径200 km内；至0.8 ms时，回击引起的电磁场几乎都在击穿阈值之下。虽然电场击穿持续时间短暂，但电子密度因雪崩电离而显著增加，图6（b）可看出电子密度增加区域分为2处。一处位于雷击通道上空80 km处，此处为准静态场引起的电子密度增强；另一处在海拔高度80～90 km间，随着电磁脉冲辐射场向外扩散至距离200 km 内。电子密度增强区域下存在微弱的电子密度减弱区域，这是因为未达到空气击穿阈值时电场对电子的加热表现为吸附作用。图6中可看出，在0.8 ms，电场减弱至击穿阈值下，但电子密度扰动区域变化不大，这可从图5 中的曲线变化看出；E/Ek比值为1.4 时，电子密度增强速率|vi-va|的值为4.67 × 108，而当E/Ek比值为0.6 时，电子密度减少速率|vi-va|的值为1.62 × 107，两者相差28.8倍，电子密度减少至常值需要的时间远比增强所需时间长。

图6 夜间击穿电场（a）及电子密度（b）随时间分布变化

3 结论

本文从麦克斯韦方程组出发，采用柱坐标二维时域有限差分方法，计入电场对电离层加热的非线性效应，基于国际标准经验模式IRI-2020、NRLMSISE-00 提供的参数和以此计算的背景电导率模式，建立了一个地闪回击电磁场在电离层传播的数值模型，并以此分析电离层电参数的昼夜差异对地闪回击电磁场传播的影响。数值模拟了电离层地闪回击产生的电磁场在电离层传播事件，并以此分析了电离层电参数的昼夜差异对传播产生的影响，主要结果如下。

1）以220 kA 峰值电流、上升特征时间为10 μs、通道高度10 km 的强地闪回击电流为源，以海南岛北部夏季白天及晚上的背景参数模拟了昼夜电离层对地闪回击强电磁场的响应。模拟结果发现：白天辐射场值在海拔60 km 处就开始衰减，而夜间辐射场至海拔80 km 处才开始衰减，夜间电离层反射高度比白天更高，回击产生的强电磁波在晚上能穿透到更高的高度。

2）空气击穿电场阈值Ek与局地空气分子数密度有关，海拔高度越高电场阈值越低。夜间地闪回击产生的电磁场相较白天更容易渗透到海拔高的区域，因此夜间高层比白天更容易产生空气击穿，产生电子密度短暂增强和瞬态发光事件等现象。