1 引言

闪电，是发生在强对流天气过程中的一类放电现象，其放电过程往往伴随着大电流和强电磁脉冲，在传播过程中会产生较强的电磁辐射效应，对地面的建构筑物和人类的生命财产安全等具有极大的潜在破坏性，因此揭示雷电电磁场辐射传播效应对公共安全以及社会经济稳定等方面具有重要意义。闪电的能量主要集中在低频和甚低频段，传播距离可达几百公里，由于在传播过程中受到地面不同电导率介质的影响，会导致其波形的幅值和相位发生偏移，在计算雷击点位置和估算峰值电流时造成较大误差。因此如何准确地评估和量化地表不同电导率对闪电电磁脉冲传播的影响，对优化站网定位和计算放电参量等问题具有重要的科学意义。

二十世纪以来，国内外众多研究学者对雷电电磁场沿地表传播特性进行了深入研究。Cooray^[1]首次提出使用有限地表阻抗表达式的方法来计算雷电回击产生的水平电场，但该方法只适用于计算距雷击点200 m内产生的水平电场。Rubinstrin^[2]在Cooray方法基础上，提供了一个适用于计算距闪电更远处水平电场的近似算法(Cooray-Rubinstein算法，简称C-R算法)。Shoory等^[3]提出了在频域中计算雷电电磁场的天线理论方法，并得出大地电导率的变化会对水平电场产生较大影响，电导率越大，水平电场越小。随着计算机技术的不断发展，时域有限差分法(FDTD)因其可计算时间和空间中任意一点的电场和磁场，而不需要近似项，近来被广泛应用于雷电电磁场模拟计算中。Baba等^[4]利用3D-FDTD研究得出平坦有损地表上建筑物存在与否对磁场影响不大，但对垂直电场影响显著。Baba等^[5]利用2D-FDTD计算闪电击中160 m和553 m高的建筑所产生的垂直电场、水平电场和磁场，结果表明相比于水平电场和磁场，雷击物存在与否对垂直电场影响较大。Zhang等^[6-7]通过2D-FDTD计算距闪电通道100~200 m处土壤垂直分层情况下的水平电场，和距闪电通道60~200 m处土壤水平分层情况下的水平电场，得出土壤分层对水平电场有一定影响。Khosravi-Farsani等^[8]用3D-FDTD计算距闪电通道100 m~2 km处水平电场，并与用C-R算法和Norton近似公式计算得出的结果进行比较。Li等^[9]建立3D-FDTD模型研究粗糙地表对雷电电磁场传播影响，模拟计算出距离闪电通道100 m内的电场。Aoki等^[10]采用FDTD研究闪电电磁脉冲在光滑有损地表传播5~200 km范围内，回击速度、电流上升沿时间和土壤电导率等参量对水平电场、垂直电场和磁场的影响。Yu等^[11]研究了闪电通道长度h、回击速度v、地面相对介电常数ε和电导率σ对水平电场的影响，结果表明与σ和v相比较之下，ε和h对水平电场的影响可忽略不计。

以上研究学者在考虑建筑物、不规则地表、土壤分层和土壤相关电气参数等因素对闪电电磁脉冲传播影响的研究中取得了令人可喜的进展。但由于闪电活动过程中常伴随着降雨，影响当地的土壤湿度，这引发了众多研究学者对土壤湿度影响闪电电磁脉冲传播方面的研究。Scott^[12]在频率10²~10⁶ Hz范围之间，对水分含量不同的样品土壤的电导率和介电常数进行测量，建立了不同湿度条件下，土壤电导率和介电常数与频率之间的经验关系参量。Longmire等^[13]基于Scott的经验关系参量，建立了土壤电参数等效模型(L-S模型)，并且提出在不同土壤湿度下，大地电导率和相对电容率是随频率变化而变化。欧阳双等^[14]和张源源等^[15]先后利用L-S模型，综合考虑土壤湿度和地表起伏，分析土壤湿度对不规则起伏地表上地闪回击电磁场传播的影响，研究得出在一定范围内土壤湿度对雷电电磁场传播有较大影响。Li等^[16]利用改进的L-S模型，研究土壤水平分层、土壤含水量和频变相关参数对闪电首次回击和继后回击在近、中和远场处产生的水平电场、垂直电场和磁场的影响。

由此可见，在影响闪电电磁脉冲传播的参数中，土壤的电特性是重要的参数之一，但之前的学者主要基于土壤湿度与电导率的经验关系进行研究讨论。本文首次将工程勘察领域研究得出的不同土壤类型及湿度和电导率计算关系式应用到雷电电磁场的模拟计算中，建立2D-FDTD模型计算水平电场、垂直电场和磁场，可更准确地研究不同土壤类型及湿度对雷电电磁场沿地表传播的影响。结果表明，在粘土、粉土和砂土这三种具有代表性的土壤类型中，雷电电磁场的水平电场分量沿粘土地表传播过程中更容易受土壤湿度影响，特别是在土壤湿度较低时，轻微的湿度变化就会对其产生显著的影响，因此土壤类型及湿度的差异对雷电电磁场沿地表传播的影响不容忽视。

2 电磁场计算方法 2.1 电流模型的选取

目前常用的雷电电流模型有双指数函数和Heidler函数模型两种，其中双指数函数模型表达式为：

$ i(0, t) = {I_0}[\exp ( - \alpha t) - \exp ( - \beta t)] $

(1)

式中I₀为通道底部电流峰值，α和β为时间常数。由于双指数函数在t=0时没有连续的一阶导数，这与实际观测到的闪电回击电流波形不符。因此，国际电工委员会(IEC)推荐使用的雷电流函数为Heidler函数，其数学表达式为：

$ i(0, t) = \frac{{{i_{01}}}}{{{\eta _1}}}\frac{{{{\left( {t/{\tau _{11}}} \right)}^2}}}{{\left[ {{{\left( {t/{\tau _{11}}} \right)}^2} + 1} \right]}}{{\rm{e}}^{ - t/{\tau _{12}}}} + \frac{{{i_{02}}}}{{{\eta _2}}}\frac{{{{\left( {t/{\tau _{21}}} \right)}^2}}}{{\left[ {{{\left( {t/{\tau _{21}}} \right)}^2} + 1} \right]}}{{\rm{e}}^{ - t/{\tau _{22}}}} $

(2)

$ {\eta _1} = \exp \left[ { - \left( {\frac{{{\tau _{11}}}}{{{\tau _{12}}}}} \right){{\left( {2 \cdot \frac{{{\tau _{12}}}}{{{\tau _{11}}}}} \right)}^{1/2}}} \right], {\eta _2} = \exp \left[ { - \left( {\frac{{{\tau _{21}}}}{{{\tau _{22}}}}} \right){{\left( {2 \cdot \frac{{{\tau _{22}}}}{{{\tau _{21}}}}} \right)}^{1/2}}} \right] $

(3)

式中i₀₁、i₀₂是回击电流幅值，τ₁₁、τ₂₁为雷电流波形的上升沿时间常数，τ₁₂、τ₂₂为雷电流波形的下降沿时间常数，Rachidi等^[17]对这些参数的取值如表 1所示，电流波形如图 1所示。

Heidler函数在t=0时对时间的导数为0，这与实际观测到的首次回击电流波形一致^[18]，此外，模型中雷电流的各个特征值都能明显地在表达式和波形图中表现出来。

2.2 回击模型的选取

雷电流回击模型是对闪电通道中电流和基电流之间的关系进行数学描述的回击规律模型，现有的回击模型分为物理模型、分布电路模型、电磁模型和工程模型这四种类型，其中工程模型运用广泛^[19]。在综合考虑计算简便性和重现电磁场特性两个方面，工程模型中的MTLL(the modified transmissionline model with linear current dacay with height)模型是研究雷电电磁场的最优模型^[20]。本文采用MTLL模型，其特点是闪电通道中的电流随通道高度增加呈线性衰减趋势，但波形形状不变；并且通道顶端回击电流为零，这样计算电磁场就忽略电流不连续性的影响。其表达式为：

$ i\left( {{z^\prime }, t} \right) = \left( {1 - {z^\prime }/H} \right)i\left( {0, t - {z^\prime }/v} \right), t \ge {z^\prime }/v $

(4)

其中，i(z', t)为回击电流沿回击通道的分布。

2.3 计算模型

本文采用的雷电电磁场计算模型是时域有限差分算法，其主要原理是对Maxwell方程在时间和空间上进行离散，即把电场和磁场在时间和空间中交错排列，将安培定律和法拉第定律中所有导数用差分形式表示，然后推导差分方程组获得递推方程组，最后递推方程组用已知场表示未知场，即计算第一个步骤的电场和磁场，使它们变为已知场，再求下一步未知的电场和磁场使之成为已知场，以此类推，直到获得在理想时间段的电磁场。

下文公式推导中Δt代表时间步长，Δr、Δz分别代表在r轴和z轴的网格时空步长，σ₀代表电导率，ε₀代表电容率，μ₀代表磁导率。在二维柱坐标系中，Yee元细胞中雷电电磁场三个分量H_φ、E_r、E_z在时间和空间上交替取样的间隔为Δt/2(表 2)。

所以水平电场E_r表示为E_rⁿ(i+1/2, j)，垂直电场E_z表示为E_zⁿ(i, j+1/2)，磁场H_φ表示为H_φ^n+1/2(i+1/2, j+1/2)，其中上标n表示第n个时间步长。由波动方程组可得：

$ - {\mu _0}\frac{{\partial {H_\varphi }}}{{\partial t}} = \nabla \times E = \left| {\begin{array}{*{20}{c}} {{{\hat a}_r}}&{{{\hat a}_\theta }}&{{{\hat a}_z}}\\ {\frac{\partial }{{\partial r}}}&0&{\frac{\partial }{{\partial z}}}\\ {{E_r}}&0&{{E_z}} \end{array}} \right| = \frac{{\partial {E_r}}}{{\partial z}} - \frac{{\partial {E_z}}}{{\partial r}} $

(5)

$ {\varepsilon _0}\frac{{\partial E}}{{\partial t}} = \nabla \times H = \frac{1}{r}\left| {\begin{array}{*{20}{c}} {{{\hat a}_r}}&{{{\hat a}_\theta }}&{{{\hat a}_z}}\\ {\frac{\partial }{{\partial r}}}&0&{\frac{\partial }{{\partial z}}}\\ 0&{r{H_\varphi }}&0 \end{array}} \right| = - \frac{{\partial {H_\varphi }}}{{\partial z}} + \frac{1}{r}\frac{{\partial {\rm{r}}{H_\varphi }}}{{\partial r}} $

(6)

整理式(5)和(6)可得：

$ {\mu _0}\frac{{\partial {H_\varphi }}}{{\partial t}} = \frac{{\partial {E_z}}}{{\partial r}} - \frac{{\partial {E_r}}}{{\partial z}} $

(7)

$ - \frac{{\partial {H_\varphi }}}{{\partial z}} = {\varepsilon _0}\frac{{\partial {E_r}}}{{\partial t}} + {\sigma _0}{E_r} $

(8)

$ \frac{1}{r}\frac{{\partial r{H_\varphi }}}{{\partial r}} = {\varepsilon _0}\frac{{\partial {E_z}}}{{\partial t}} + {\sigma _0}{E_z} $

(9)

对式(7)差分并整理得：

$ \begin{array}{c} H_\varphi ^{n + 1/2}(i + 1/2, j + 1/2) = {H_\varphi }^{n - 1/2}(i + 1/2, j + 1/2) + \frac{{\Delta t}}{{{\mu _0}\Delta r}}\left[ {{E_z}^n(i + 1, j + 1/2) - E_r^n(i, j + 1/2)} \right]\\ - \frac{{\Delta t}}{{{\mu _0}\Delta z}}\left[ {E_r^n(i + 1/2, j + 1) - E_r^n(i + 1/2, j)} \right] \end{array} $

(10)

对式(8)差分并整理得：

$ {E_r}^{n + 1}(i + 1/2, j) = A \cdot E_r^n(i + 1/2, j) - \frac{B}{{\Delta z}} \cdot H_\varphi ^{n + 1/2}(i + 1/2, j + 1/2) - \frac{B}{{\Delta z}} \cdot H_\varphi ^{n + 1/2}(i + 1/2, j - 1/2) $

(11)

对式(9)差分并整理得：

$ \begin{array}{c} {E_z}^{n + 1}(i, j + 1/2) = A \cdot {E_z}^n(i, j + 1/2) + \frac{B}{{i\Delta r}} \cdot (i + 1/2)H_\varphi ^{n + 1/2}(i + 1/2, j + 1/2)\\ - \frac{B}{{i\Delta r}} \cdot (i - 1/2)H_\varphi ^{n + 1/2}(i + 1/2, j + 1/2) \end{array} $

(12)

其中：$ A = \frac{{2{\varepsilon _0} - {\mathit{\sigma }_0}\Delta t}}{{2{\varepsilon _0} + {\mathit{\sigma }_0}\Delta t}}, B = \frac{{2\Delta t}}{{2{\varepsilon _0} + {\mathit{\sigma }_0}\Delta \mathit{t}}}$

2.4 电磁场传播模型

图 2为雷电电磁场沿地表传播的计算模型，ε₀和ε₁分别代表空气和大地的电容率，μ₀和μ₁分别代表空气和大地的磁导率，σ₁代表大地的电导率，其计算空间为5 000 m×5 000 m的一柱切面，计算时间步长Δt=1.66 ns，吸收边界采用Mur一阶边界条件^[21]。假设大地为均匀有损土壤，闪电通道高度设为H=8 km，观测点的高度设为h=10 m，雷电回击速度v设为光速c的一半，即v=c 2=1.5×10⁸m/s^[11]。d为观测点与闪电通道距离，为探究不同距离处雷电电磁场的变化，本文取4个观测点，分别为50 m、100 m、200 m和1 km。

2.5 土壤电气参数变化模型

孙旭等^[22]基于多孔介质模型分析土壤的导电机制，根据控制变量的研究思想，采用强迫电流法，通过构建土壤导电机制模型研究土壤电导率与湿度的关系。为了使测量结果具有普遍性，他们根据国家标准《土的分类标准》采集了三种具有代表性的土壤，分别为粘土、粉土和砂土，并将土壤烘干击碎磨细后，用孔径0.5 mm的筛子过筛，选取颗粒较小的土样按照式(13)进行不同湿度土样的配置：

$ {m_w} = \frac{m}{{1 + 0.01{w_0}}} \times 0.01\left( {{w^\prime } - {w_0}} \right) $

(13)

式中m_w为土样所需的加水量，m为土样质量，w₀为风干含水率，w'为制备土样要求的含水率。为保证测量数据的准确性，每种湿度土样备制3个模型，每个模型测试5次以上，用最小二乘法对实验数据进行处理，从图 3可以看出，在土壤湿度低于30%的情况下的拟合效果显著。

孙旭等^[22]鉴于以上分析，对低含水率的电导率数据单独进行线性拟合，得到粘土、粉土和砂土的土壤湿度与电导率之间的线性拟合公式如下：

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _{{\rm{caly\;粘土}}}} = 0.005\;727p - 0.025\;120}\\ {{\sigma _{{\rm{silt\;粉土}}}} = 0.003\;111p + 0.006\;215}\\ {{\sigma _{{\rm{sand\;砂土}}}} = 0.001\;437p - 0.002\;342} \end{array}} \right. $

(14)

上式中p为土壤湿度，σ表示土壤电导率，下标表示不同类型的土壤。

3 数值计算结果 3.1 不同土壤类型下的水平电场

我国地区土壤湿度主要集中在5%~25%之间^[23]，通过2.5节中土壤湿度与电导率的计算公式，表 3给出了粘土、粉土和砂土在各湿度下所对应的电导率值。为研究水平电场沿这三种类型土壤地表传播的特性，图 4和图 5给出了观测点在200 m处，土壤湿度分别为5%、6%、10%、15%和25%时水平电场变化趋势图。

从图 4中可以看出，这三种土壤的水平电场强度幅值随着土壤湿度的增大而减小，因为湿度越大，电导率越大，电场的高频分量在增加^[24]。相比于粉土和砂土，粘土时的水平电场受土壤湿度影响更为显著，其强度幅值在土壤湿度从5%到25%的变化过程中波动最明显，幅值衰减最大。在土壤湿度较低时，轻微的湿度变化就会对粘土的水平电场产生显著的影响。从图 5中可以看出，土壤类型为粘土时，水平电场强度幅值在土壤湿度从5%到6%变化过程中跃变到小于砂土的水平电场强度幅值，当湿度增大到15%时跃变到小于粉土的水平电场强度幅值。由此可见，土壤类型的不同对水平电场的影响较为显著，因为土壤的电特性主要取决于土壤湿度、土壤颗粒表面的吸附特性及其相互连结作用，其中土壤湿度是影响土壤电导率的最主要因素。土壤属于多孔介质，其湿度取决于当地水文条件，以及土壤孔隙溶液中离子含量与矿物质在土壤中溶解度和电离特性有关^[25]。粘土含沙量比砂土和粉土少，其渗水速度慢，保水性的能力好，并且粘土中土壤与土壤之间孔隙比较多，当发生降雨时，土壤中水分增多，溶液会慢慢地填满土壤中各个孔隙，形成电通路，土壤的导电能力增强，电导率变大，电场的高频分量就会增加^[24]。

表 4、表 5和表 6给出了这三种土壤在各湿度变化区间内土壤湿度每增加百分之一，水平电场峰值减小的幅度，其中Δp表示各土壤湿度区间的差值，Δ表示这4个湿度区间中土壤湿度每增加百分之一水平电场峰值所减小的幅度。从表格中Δp和Δ可以看出，在土壤湿度较低时，轻微的湿度变化就会对水平电场强度幅值产生显著影响，特别是在土壤类型为粘土时，土壤湿度从5%到6%变化过程中，水平电场峰值的波动幅度达到39.86%。随着土壤湿度的增大，粘土、粉土和砂土的Δ值越来越小，这是因为当土壤水分持续增加时，土壤中可以填满的孔隙逐渐达到饱和，随之形成的导电通路也逐渐饱和，土壤电导率的上升趋势减缓^[26]，此时土壤湿度的改变就对水平电场强度幅值的影响减弱了。

3.2 土壤湿度对各场分量的影响

由3.1节分析得出土壤类型为粘土时，水平电场受土壤湿度影响较为显著，所以本节计算分析土壤类型为粘土时，土壤湿度对各观测点处雷电电磁场磁场分量H_ⱷ，垂直电场分量E_z和水平电场分量E_r的影响。从图 6、图 7和图 8中可以看出，在同一个观测点，土壤湿度的改变对水平电场有显著的影响，其幅值随着土壤湿度的增加而显著减小，对磁场和垂直电场的影响基本可以忽略，磁场和垂直电场强度幅值随时间变化的趋势几乎相同^[10]。随着观测距离的增加，磁场、垂直电场和水平电场的强度幅值显著减小，波头上升沿时间显著变慢。这是因为电磁场在传播过程中，低频分量不容易衰减可以传播到很远的距离，而高频分量衰减迅速，所以时域场峰值减小、波形上升沿时间增大。

为得出土壤湿度的改变对磁场和垂直电场影响的程度等级，量化磁场和垂直电场峰值数据(表 7)，其中ΔE_z1和ΔH₁分别代表土壤湿度为5%和15%时垂直电场峰值的差值和磁场峰值的差值，ΔE_z2和ΔH₂分别代表土壤湿度为15%和25%时垂直电场峰值的差值和磁场峰值的差值。从表中数据可以看出，ΔE_z1和ΔE_z2随土壤湿度的增加而减小，ΔH₁和ΔH₂值却近似相等，所以相比于磁场，垂直电场更容易受土壤湿度改变的影响。

为更直观地看出土壤湿度对水平电场波头上升沿时间的影响，表 8列出了水平电场在土壤湿度为5%、15%和25%，观测点为50 m、100 m、200 m和1 km处，场峰值所对应的上升沿时间。对比表中数据可以得出，随着土壤湿度的增大，水平电场波头上升沿时间变慢。不同土壤湿度对雷电电磁场的水平电场分量有显著的影响，对垂直电场和磁场分量影响较小^[10]。这是因为土壤电导率随着土壤湿度的改变而改变，由于土壤电导率对雷电电磁脉冲高频信号有衰减作用，不同的电导率对雷电电磁脉冲衰减作用不同，导致雷电电磁场各分量幅值的差异^[26]。雷击电磁场主要分为静电场、感应场和辐射场，垂直电场在近场主要表现为静电场，辐射场很小，磁场在近场主要表现为感应场，而土壤电导率对辐射场有所衰减，所以土壤湿度的改变对近场处的垂直电场和磁场影响非常小。

4 讨论和总结

不同土壤类型中土壤颗粒组成比例不同，其渗水和保水能力差异等因素会导致土壤的电特性不同，因此土壤是影响闪电电磁脉冲沿地表传播特性的重要参数之一^[16]。闪电发生时常伴随着降雨，这使土壤电气特性有所改变，其中土壤电导率的改变对雷电电磁场影响尤为显著^[9]。闪电定位系统利用到达时间差法(Time of arrival algorithm，简称TOA)来计算闪电的位置，其原理是根据闪电辐射的电磁波到达系统中各个测站的绝对时间差来确定闪电的位置。然而，闪电电磁脉冲在空间传播过程中往往受到大气折射、地形和土壤电气特性等因素影响，造成波形相位和幅值的偏移，这会对计算雷击点位置和估算峰值电流造成较大误差。

本文首次将工程勘察领域研究得出的土壤湿度和电导率计算关系式应用到雷电电磁场模拟中，基于Heidler雷电通道基电流函数模型和MTLL回击模型，在Mur一阶边界条件下，利用2D-FDTD模拟粘土、粉土和砂土这三种土壤在不同湿度条件下的水平电场，并计算土壤类型为粘土时，观测点在50 m、100 m、200 m和1 km处，土壤湿度为5%、15%和25%条件下的水平电场E_r、垂直电场E_z和磁场H_ⱷ。

(1) 在粘土、粉土和砂土这三种土壤中，粘土的水平电场受土壤湿度影响最为明显，即土壤湿度增大，场强度幅值减小，波头上升沿时间变慢。特别是土壤湿度较低时，轻微的湿度变化就会对其强度幅值产生大幅度波动，即土壤湿度从5%到6%变化过程中，水平电场峰值的波动幅度达到39.86%。

(2) 当土壤水分增加到一定值时，土壤电导率的上升趋势减缓，此时土壤湿度的改变对水平电场的影响减弱。

(3) 土壤湿度对雷电电磁场三个分量影响的程度等级为E_r > E_z > H_φ，即水平电场受土壤湿度影响显著，垂直电场和磁场受土壤湿度的影响可基本忽略。

(4) 土壤湿度不变时，随着观测点与闪电通道距离的增大，沿地表传播的电磁场强度幅值显著减小，波头上升沿时间显著变慢。

本文研究结果表明，土壤类型和湿度不同，一定程度上会影响闪电电磁脉冲沿地表的传播特性，主要表现为场幅值衰减和波头上升沿时间变慢两个方面，而TOA闪电定位是对接收到的闪电电磁脉冲参量进行放电参数的反演和闪电定位的估算，所以在对放电参数进行反演计算时，应考虑电磁场峰值衰减引起的误差，在估算闪电始发位置时，要剔除因各地区之间土壤类型的差异和土壤湿度的改变引起的上升沿时间误差，这对优化闪电站网定位和计算放电参量具有重要的科学意义。未来，我们将进一步探讨一定倾角入射的云闪脉冲衰减特征与土壤湿度之间的关系。