引言

雷电灾害对人类社会构成巨大威胁，被联合国减灾十年委员会列为“十大最严重的自然灾害之一”。雷电灾害在中国的总体分布是东南多、西北少, 南方地区雷电灾害占总雷灾的82.98%，雷灾伤亡人数也占到总雷灾伤亡人数的82.94%(Holle，2016；殷启元等，2019；Yin et al，2021)。随着中国经济和社会发展的进步，国内不同区域和不同行业都迫切需要对雷电进行监测和预报，以保护重要设施、区域和活动(许小峰，2004；吕伟涛等，2009；马颖等，2009；刘维成等，2015；张烨方等，2021)。大气电场强度是大气物理和大气电学的基本参数(Odzimek et al，2018；Bernard et al，2020；郄秀书等，2024)，雷电活动常伴随着大气电场的变化，这些变化可以用来进行雷电风险预警(孟青等，2005；方俏娴等，2023)。近年来基于大气电场的雷电风险预警研究发展迅速，预警方法大致分为时域方法和频域方法。

时域方法通常采用易于实现的电场阈值法进行雷电风险预警，Murphy et al(2008)对大气电场数据采用阈值判定法进行了初步雷电风险预警的探索；Aranguren et al(2009)对比分析了阈值法和极性反转法，发现监测电场数据发生首次极性反转更能有效地预报闪电的发生，预报的准确率达到47%。Srivastava et al(2015)结合大气电场与马尔可夫模型，预警准确率达到了66.4%，但其误报率也达到了59.7%，故无法投入实际应用。基于时域分析的方法准确率低，且忽略了电场信号的固有物理特性，未能深入探索电场信号特征与雷电之间的关系(Montanya et al，2004)。

频域方法相比较时域方法预警效果虽得到了一定提升，但对于大气电场资料的挖掘仍旧不够充分(李颖等，2013；徐栋璞等，2013)。康海龙等(2014)通过短时傅里叶变换分析发现能量特征可作为雷电风险预警的指示特征；Lu et al(2016)利用希尔伯特-黄变换分析发现雷电生成过程中其信号高频能量逐渐升高，从而提出使用信号的高频能量占比进行雷电风险预警；刘刈(2021)利用希尔伯特-黄变换分析重庆北碚站不同天气过程中大气电场资料的震荡，并开展雷电风险预警方法研究。但是，大气电场资料采集过程中含有空间噪声信号，影响其真实反映天气现象，上述研究均未考虑噪声对预警效果的影响(Yang et al，2021)。

大气电场数据作为一种非平稳信号，其安装位置以及周围环境存在很多不利因素影响着真实大气电场特征的识别。小波分解能够对输入信号的细节分量进行快速分类和分析，庞金凤等(2019)通过小波分解方法分析2008—2016年沙尘天气发生时序，分解为平稳性和非线性。数字形态滤波器是一种重要的非线性滤波技术，对一维信号处理能够有效地抑制噪声的干扰，对混噪中有用信息质量的增强效果显著。刘金江等(2007)利用形态学和小波包理论相结合的方法用于心电图检测方面的判定。李晓梅等(2018)利用形态学开-闭重构方法在背景模型中抑制噪声来进行精准分割，精度可以达到90%。

在本文中，为利用好大气电场信号特征且提高雷电风险预警精准度，基于增强经验小波变换和自适应平滑滤波法(enhanced empirical wavelet transform+adaptive Savitzky-Golay, EEWT-ASG)，提出了一种结合时频域特征和一维形态学的雷电风险预警方法。该方法利用EEWT-ASG对大气电场数据进行去噪平滑处理和时频域特征提取，再通过数字形态学进行全局趋势分析计算，最终实现精准雷电风险预警目的。该方法相比较单独基于时域或频域特征进行预警的方法，提升了准确率，同时也降低了漏报率和误报率。

1 数据与方法 1.1 数据

本文大气电场数据(atmosphere electric field, AEF)来源于安装在广州塔平台(500 m高度)的大气电场仪，该设备采用高灵敏、低功耗的MEMS电场传感器敏感芯片技术(Yang et al，2015)。经过试验比对，该款大气电场数据与传统的场磨式电场仪探测结果相似，能够满足对雷暴天气的实时监测(Zhou et al，2016)。

本文参考“两区域法”(Murphy et al，2008；Xu et al，2022)，即关注区(area of concern, AOC)半径10 km、警告区(warning area, WA)半径20 km进行选取样本和效果验证。如图 1，AOC是目标点安全保护范围，WA是AOC外的一个警戒区域。选取了2021年和2022年距离大气电场仪20 km内83天有雷暴过程和123天非雷暴过程对应的大气电场仪全天(00:00—24:00)(北京时，下同)资料。本文对雷暴的判别依据是指该范围内是否有雷电发生，包括云闪和地闪。

1.2 方法

在本文中，提出基于EEWT-ASG开展对大气电场信号的噪声处理和基于数字形态学的全局趋势计算方法，提升雷电风险预警方法的鲁棒性和普适性。

1.2.1 EEWT-ASG

Yang et al(2015)和Lu et al(2016)指出去噪是研究AEF的首要任务，并提出基于迭代思想的信号分解方法处理大气电场噪声。但由于该方法处理复杂信号用时较长，并不适用于雷电这种实时性极高的风险预警任务。Gilles(2013)提出经验小波变换(empirical wavelet transform, EWT)，其相比较经验模态分解的运算速度更快，然而EWT在合理处理频谱分割时会造成过度分割(Zhang et al，2023)，导致增加不必要的冗余信息。Hu et al(2017)对频谱分割方法进行了研究提升，在EWT的基础上提出了一种增强型的经验小波变换方法(EEWT)。

本文使用EEWT方法进行AEF信号的分解，AEF信号及分解后的k个模态信号分别记为D(n)和D_i(n)，i=1, 2, …, k。

而后，采用自相关分析法来识别D_i(n)中包含噪声较多的模态信号。本文将理想高斯白噪声和D_i(n)的自相关函数计算结果分别记为R_N和R_Di。

使用皮尔逊算子(P)计算R_Di与R_N的相似性(Hu et al, 2017)：

$ P=\frac{\operatorname{COV}\left(R_{D_i}, R_N\right)}{\sigma_{R_{D_i}} \sigma_{R_N}} \quad i=1, 2, \cdots, k $

(1)

式中，COV、σ分别代表协方差和标准差。P越接近1，表示R_Di与R_N之间的相关性越好，则对应D_i(n)中包含噪声的可能性就越大。本文定义包含噪声的模态信号D_N(n)如下：

$ D_N(n)= \begin{cases}D_i(n) & P_i>0.75 \\ 0 & P_i \leqslant 0.75\end{cases} $

(2)

Savitzky-Golay滤波器设计过程中有两个关键参数，即数据窗口大小和多项式阶数。其中，数据窗口如果过大，会导致信号中部分有用信息的丢失，如果过小则影响过滤信号效果；而多项式阶数的选择，如果过大可能会引入新的噪声，过小会导致信号过度平滑而失真(Lu et al，2016；Muthusivagami，2020)。由于大气电场信号为非平稳信号，在信号平滑过程中采用固定的数据窗口大小和多项式阶数容易造成信号中原有信息的丢失。

因此，本文根据Yang et al(2021)研究，多项式拟合阶数选为3阶，同时利用广义无偏估计均方差(generalized unbiased estimate of MSE，GUE-MSE) 实现固定阶数动态调整数据窗口大小(John et al，2021)，对式(2)中包含噪声的模态信号D_N(n)进行过滤。滤波结果如图 2所示。图 2a~2d分别为原始采集波形、基于GUE-MSE方案动态调整的数据窗口波形、经过窗口滤波后的波形、被滤波器滤掉的噪声信号波形。由图可见，该方法在去噪的同时，能够较为真实地还原大气电场信号。

1.2.2 基于形态学的全局趋势计算方法

受趋势计算(Gilles，2013；Hu et al，2017)的启发，为准确描述一段时间内信号特征，在数字形态学(Morpho)的基础上提出以下计算方法。其中一维膨胀和腐蚀计算方法就是根据滑动窗口取数据序列的最大和最小值，滑动窗口的长度是通过计算数据局部极大值之间的欧几里得距离的最大值决定的，Morpho是计算闭运算和开运算的算术平均值。

在进行膨胀和腐蚀计算前，为了保证处理前后数据大小一致，需使用镜像扩展法(Gilles，2013)对数据进行预处理。Morpho的计算结果如图 3中黑色曲线所示，可以看到该曲线虽然可以显示该段信号的趋势变化，但在计算结果中存在与时间轴平行的复杂区间，本文定义该区间为平顶区间。为了更好描述平顶区间的变化趋势，本文对Morpho的计算结果进行修正。

首先计算变化率ΔMorpho(n)，在此基础上定义局部趋势函数(Local)：

$ \text { Local }= \begin{cases}\mathrm{Upw} & \Delta \operatorname{Morpho}(n)>0 \\ \text { FT } & \Delta \operatorname{Morpho}(n)=0 \\ \text { Downw } & \Delta \operatorname{Morpho}(n)<0\end{cases} $

(3)

式中：Upw、Downw、FT分别指上升区间、下降区间、平顶区间。

而后，对平顶区间进行进一步细分：

① 如果其左、右边均为上升区间，则该平顶区间标记为上升区间；

② 如果其左、右边均为下降区间，则该平顶区间标记为下降区间；

③ 其他情况，记为复杂区间。

基于局部趋势函数的结果，本文定义基于全局趋势(Global)的判断规则[式(4)]，进一步修正对变化趋势的判断。

$ \begin{gathered} \text { Global }= \\ \left\{\begin{array}{lc} \text { Uptr } & \frac{L(\text { Upw })}{L(\text { Upw })+L(\text { Downw })}>0.60 \\ \text { Downtr } & \frac{L(\text { Upw })}{L(\text { Upw })+L(\text { Downw })}>0.60 \\ \text { Unc } & \text { 其他 } \end{array}\right. \end{gathered} $

(4)

式中：Uptr、Downtr分别为修正后的上升区间、下降区间，Unc为复杂区间；L指上升(下降)区间的相应时长。

经Morpho修正后的结果如图 3所示。该判断规则相比较使用斜率判断全局趋势而言，分析更为细致和精准。同时，由于该方法滑动窗口的大小取决于数据局部极大值之间欧几里得距离的最大值，因此，在保证对毛刺信号有一定过滤效果的同时，可保障不会由于局部信号的大幅度抖动而对全局趋势的判定造成影响。

2 雷电风险预警方法 2.1 时频谱特征统计

雷暴过程的大气电场变化及时频谱如图 4所示。其中，红色及绿色虚线分别表示结合闪电定位和雷达资料，经人工研判确定的预警时间及解除预警时间。由图可见，在预警时间附近(红色虚线)大气电场信号发生剧烈抖动，同时谱宽显著增大。因此统计谱宽B(n)和能量差Diff的变化特征，以提取雷电风险预警判定标准。

首先分析谱宽的变化情况。对B(n)计算一阶后向差分，记为ΔB(n)。AEF数据的计算长度为10 min，统计了雷暴和非雷暴过程ΔB(n)的平均值、最大值、最小值及其标准差，分别记为|ΔB_avg|、|ΔB_max|、|ΔB_min|、STD，结果如表 1所示。其中|ΔB_min|均为0，故表中未逐项列出。本文将雷暴过程按照获取的预警信息从时间角度分为四个阶段：预警前、预警后、解除预警前和解除预警后。从表 1可以发现，|ΔB_avg|和STD可以明显区分雷暴和非雷暴过程。因此，选择这两个指标作为雷电风险预警特征是比较合理的。

进一步对大气电场能量差Diff的变化进行分析。设AEF信号及去除低频成分后AEF信号的能量分别为E(n)和E′(n)，触发预警和解除预警时刻分别为n_w和n_dw，E(n)的能量差Diff如式(5)。E′(n)的能量差Diff′计算方法类似。对Diff和Diff′统计其相应的均值、最大值和最小值，分别记为Diff_avg(Diff′_avg)、Diff_max(Diff′_max)和Diff_min(Diff′_min), 详见表 2。

$ \operatorname{Diff}=\left[\begin{array}{ccc} E\left(n_{\mathrm{w}}-300, n_{\mathrm{dw}}-300\right) & \cdots & E\left(n_{\mathrm{w}}-300, n_{\mathrm{dw}}+300\right) \\ E\left(n_{\mathrm{w}}-299, n_{\mathrm{dw}}-300\right) & \cdots & E\left(n_{\mathrm{w}}-299, n_{\mathrm{dw}}+300\right) \\ \vdots & & \vdots \\ E\left(n_{\mathrm{w}}+300, n_{\mathrm{dw}}-300\right) & \cdots & E\left(n_{\mathrm{w}}+300, n_{\mathrm{dw}}+300\right) \end{array}\right] $

(5)

结合表 2，可以发现n_w和n_dw之间的能量差总体来说很小。同时，考虑到在非雷暴过程中AEF信号存在缓慢的上下波动(即低频信号)，在雷暴预警中本文将其视为干扰。在去除低频信号能量后，Diff′_min仅为1 dB左右，所以本文选择E′(n)作为主要的解除预警特征。对比原始AEF信号和滤波后的AEF信号，可以看到在经过滤波后Diff′_avg和Diff′_max的范围都有一定程度的缩小。这说明滤波后的AEF信号在能量上表现出n_w和n_dw之间更好的相关性，这种相关性的提升对雷电风险预警有着积极作用。

2.2 雷电风险预警方法

在2.1节统计分析基础上，本文的雷电风险预警方法可分为触发预警和解除预警两部分，下文中提及的“趋势”使用1.2.2节方法进行计算。同时，定义E′(n)在对应时长内的均值为E′_avg。

触发预警条件：

① 20 min范围的AEF信号B(n)和E′(n)呈全局上升趋势

② 10 min范围的E′_avg>-20 dB、STD>0.5 kV·m^-1、|ΔB_avg|>0.5 mHz

解除预警条件：

本文设计了两个方案并行判断来共同实现自动解除雷电风险预警。其中任意一个方案的条件被满足或两个方案同时满足，即可解除当前预警。两个方案分别如下：

方案1：

① 20 min范围的AEF信号的B(n)和E′(n)呈全局下降趋势;

② 10 min范围的|E′_avg-E′_w| < 1 dB、STD < 1 kV·m^-1。

方案2：

①|ΔB_avg|连续10次判断的总和小于0.1 mHz;

② 10 min范围的STD < 1 kV·m^-1。

方案1中E′_w为预警时刻之前的10次判断的E′_avg最小值。E′_w随着每次判断为预警时E′_avg的不同而不同。预警方法流程如图 5所示。本文考虑到实际的业务应用需要：①使用EEWT-ASG对最近20 min的大气电场数据进行分解识别噪声分量，并进行去噪处理重构数据；②利用小波变换提取重构数据的时域和频域特征，通过一维形态学获取对应特征的趋势；③根据上述的特征状态，判断是否需要发出预警或者取消预警；④输出当前的判断结果，等待下一次数据输入。

3 效果检验

本文的效果检验是指有效预警为在一个警告周期内AOC内首次地闪之前发出的预警；未能预警为在一个警告周期内AOC内首次地闪之后发出的预警或未触发预警；错误预警为在一个警告周期内AOC内没有出现地闪。根据相关研究(Barnes et al，2009; Meng et al，2019；Xu et al，2022)，本文采用以下5个指标来评价预警行为的性能：准确率(POD)、漏报率(MAR)、误报率(FAR)、临界成功指数(CSI)和平均预警提前量。

采用上述指标，对比本文提出的预警方法与Lu et al(2016)、Zeng et al(2013)的EFAI和EFDI方法，结果见表 3：

通过表 3可发现，对于雷暴过程中，本文提出的方法在保证了较高的POD的情况下(77.11%)，同时也具备最低的FAR(40.00%)和最高的CSI(0.51)性能，平均预警提前量也达到了22.27 min。Zeng et al(2013)所提出的EFDI方法虽然POD表现最好(91.57%)，但其FAR却达到了60.00%，这对于雷电风险预警业务使用不具备可行性。在非雷暴过程中，本文所提出的方法POD达到了90.24%，虽略逊于Lu et al(2016)的方法(98.26%)，但Lu et al(2016)的方法在雷暴过程中的POD表现却是最低的(18.07%)，也无法达到准确雷电风险预警目的(Cooper and Holle, 2019)。综合对比，本文提出的雷电风险预警方法优于其他作者的结果，具备雷电风险预警业务使用可行性。

图 6给出了利用本文方法的雷电风险预警提前量分布图，可看出预警提前量绝大部分处于0~40 min区间，其中分布在0~20 min区间的有62.5%，完全能够满足大多数雷电高影响行业需求。

目前大多数雷电风险预警方法未考虑预警解除问题，本文定义AOC内最后一次出现地闪时刻为一次雷暴过程结束的标志时间。同时，本文统计了成功提前预警的64个雷暴过程的预警解除滞后期，具体分布见图 7。从图可见，81.3%(52/64)的可以完成预警解除，且54.7%(35/64)的过程解除滞后时间处于0~40 min，平均滞后时间为32 min，这与雷电安全“双30”法则(Cooper and Holle, 2019)相近。因此，本文所提出的解除雷电风险预警方法可满足雷电高影响行业的安全风险解除要求。

4 结论与讨论

本文提出一种结合时频域特征和一维数学形态学的雷电风险预警方法，实现了面向目标点区域的准确雷电风险预警触发和解除。主要结论如下：

(1) 通过设计EEWT-ASG滤波，完成对所选雷暴和非雷暴过程的信号重构，根据时频谱特征统计，发现信号能量与雷暴过程的生消发展相关，并提出了使用时频域特征的|ΔB_avg|和STD作为雷电和非雷电过程的预警判断特征，AEF信号的能量差作为判断解除预警的特征。

(2) 通过引入形态学Morpho计算的全局趋势作为预警特征，从AEF数据的整体走势去做预警判断，区别于前人单纯阈值判断的做法，提升了方法的鲁棒性和普适性。

(3) 通过对选取样本进行效果检验，本文所提出的雷电风险预警方法在雷暴过程中POD达到了77.11%，同时也具备最低的FAR(40.00%)和最高的CSI(0.51)性能，平均预警提前量也达到了22.27 min；非雷暴过程中POD达到了90.24%；大部分预警解除滞后时间处于0~40 min，平均滞后时间为32 min。

该雷电风险预警方法仍存在一定不足，雷电风险预警的准确率还有提升空间，漏报率和误报率仍有降低空间。下一步可以考虑实现雷电高影响用户所需的分级预警能力，综合考虑真实三维高精度闪电能量信息，以及增添雷达资料和大气电场组网资料获取综合信息进行雷暴云团定位及预警工作。