1. Introduction

2. Method

2.1. EEWT-ASG

In 2021, Yang emphasized the importance of denoising for the study of AEFs and proposed a complementary ensemble empirical mode decomposition with adaptive noise and a Savitzky-Golay filter (CEEMDAN-SG) for AEF noise reduction [14]. However, CEEMDAN, being an iterative signal decomposition method, demands a considerable amount of time to process intricate signals [15], rendering it unsuitable for real-time warning requirements. In 2013, Gilles presented an empirical Wavelet transform (EWT) [16], which offers faster processing than empirical mode decomposition (EMD). Nevertheless, the challenge with EWT lies in executing spectrum segmentation. Although Gilles proposed several solutions to this issue [17], these methods still led to excessive spectrum segmentation, resulting in redundancy in the decomposition outcomes. In 2017, Hu investigated Gilles’ spectrum segmentation methods and introduced an enhanced empirical Wavelet transform [18].

In this study, we employed EEWT to carry out the decomposition of AEF signals, with the decomposition outcomes displayed in Figure 1. It is evident that the original AEF signal is partitioned into four modes, each displaying a progressively increasing frequency. We designated the AEF signal and the decomposed k component signals as $D\left(n\right)$ and $D_i\left(n\right)$, $i\in [1,k]$, respectively, which maintain the following relationship:

$$D\left(n\right)=\sum _{i=1}^k{D}_i\left(n\right).$$

(1)

Indeed, the value of k is dictated by the inherent complexity of the processed AEF signal, which stems from the adaptive spectral analysis performed by the EEWT. The more intricate the frequency components of the AEF signal, the greater the number of decomposed modes. In order to ascertain whether each component contained noise, we employed an autocorrelation analysis, inspired by Yang, to make this determination[14]. We computed the normalized autocorrelation function of $D_i\left(n\right)$ as follows:

$$R_{D_i}\left(m\right)=\frac{{\sum }_{n=-\infty }^{\infty }{D}_i\left(n\right)D_i(n+m)}{{\sum }_{n=-\infty }^{\infty }{D}_i\left(n\right)D_i\left(n\right)}.$$

(2)

At this time, the normalized autocorrelation function for all components can be obtained as :

$$R_D\approx \left[\begin{array}{ccc}{R}_{D_1}& R_{D_2}\cdots & R_{D_k}\end{array}\right].$$

(3)

Furthermore, the normalized autocorrelation function of the ideal Gaussian white noise was computed and denoted as $R_{Noise}$. Subsequently, the Pearson correlation coefficient [19] was employed to gauge the similarity between $R_{D_i}$ and $R_{Noise}$, as follows:

$$Pearson=\frac{COV(R_D,R_{Noise})}{{\sigma }_{R_D}{\sigma }_{R_{Noise}}}=\left[\begin{array}{ccc}Pearso{n}_1& Pearso{n}_2\cdots & Pearso{n}_k\end{array}\right]$$

(4)

Figure 1. The modes extracted by EEWT are displayed, with (a) representing the original signal and (b) through (e) illustrating the decomposition modes.

Figure 1. The modes extracted by EEWT are displayed, with (a) representing the original signal and (b) through (e) illustrating the decomposition modes.

Among them, $COV$ and $\sigma$ denote the covariance and standard deviation, respectively. The range of $Pearso{n}_i$ is between -1 and 1. $Pearso{n}_i$ closer to 1 means that the correlation between $R_{D_i}$ and $R_{Noise}$ is better, then the more probable that $R_{D_i}$ contains noise. In this study, the component signal $D_{Noise}\left(n\right)$ that contains noise is defined as follows:

$$D_{Noise}\left(n\right)=\left\{\begin{array}{cc}{\displaystyle D_i\left(n\right),}\hfill & \mathrm{if}Pearso{n}_i>0.75\hfill \\ {\displaystyle 0,}\hfill & \mathrm{if}Pearso{n}_i\le 0.75\hfill \end{array}\right.$$

(5)

An Adaptive Savitzky-Golay (ASG) filter is employed for $D_{Noise}\left(n\right)$ to smooth it [20]. The ASG filter addresses the issue of selecting two key parameters in the SG filter: the size of the data window and the polynomial degree. If the data window is excessively wide, it may result in the loss of valuable information in the signal. Conversely, if the data window is too narrow to effectively filter the signal, opting for a large polynomial degree might introduce new noise, while selecting a small polynomial degree could cause distortion due to oversmoothing of the signal [14,21]. Given that the AEF signal is non-stationary, fixing both the data window size and the polynomial degree in signal smoothing could lead to the loss of original information in the signal. In [14], Xu utilized an SG filter to address the noise in AEF and set the polynomial degree and data window size to 3 and 7, respectively. Building on Xu’s work, we adopted a polynomial degree of 3 and employed the G-FL scheme to dynamically adjust the data window size [20]. The denoising results are depicted in Figure 2.

Figure 2. The denoising results of EEWT-ASG are presented, with (a) representing the original signal, (b) illustrating the denoised signal, (c) depicting the dynamic adjustment of the data window size using the G-FL scheme, and (d) showcasing the noisy signal.

Figure 2. The denoising results of EEWT-ASG are presented, with (a) representing the original signal, (b) illustrating the denoised signal, (c) depicting the dynamic adjustment of the data window size using the G-FL scheme, and (d) showcasing the noisy signal.

2.2. Calculation of global trend based on Morpho

To delineate the global trend of the signal over time, drawing inspiration from Hu and Gilles’ trend computation [18,22], we put forth the following method to acquire the global trend based on Morphology.

$$Dilate\left(n\right)=MA{X}_{k\in A_s}\left(X\left(k\right)\right)$$

(6)

$$Erode\left(n\right)=MI{N}_{k\in A_s}\left(X\left(k\right)\right)$$

(7)

Where X is the data sequence and $A_s$ is a sliding window of size s. Euclidean distances between the local maxima of the data are calculated and denoted as $X_{Localmax}$.

$$s=MAX\left(X_{Localmax}\right)$$

(8)

$$Cl\left(n\right)=Erode\left(Dilate\right(n\left)\right)$$

(9)

$$Op\left(n\right)=Dilate\left(Erode\right(n\left)\right)$$

(10)

$$morpho\left(n\right)=\frac{Cl\left(n\right)+Op\left(n\right)}{2}$$

(11)

Prior to the Dilate and Erode calculations, data must undergo preprocessing using the mirror expansion method to ensure consistent data size [18]. The results of Morpho are rectified through the following steps. Firstly, the first-order backward differentiation of $morpho\left(n\right)$ is calculated and denoted as $\Delta morpho\left(n\right)$. Subsequently, the upward and downward regions in $morpho\left(n\right)$ are identified as:

$$Trend=\left\{\begin{array}{cc}{\displaystyle UpwardRegion,}\hfill & \mathrm{if}\Delta morpho>0\hfill \\ {\displaystyle FlatTop,}\hfill & \mathrm{if}\Delta morpho=0\hfill \\ {\displaystyle DownwardRegion,}\hfill & \mathrm{if}\Delta morpho<0\hfill \end{array}\right.$$

(12)

Ultimately, the flat top is rectified based on the following three criteria:

• the flat top is classified as the upward region if both the left and right sides of the flat top are within the upward region;

• the flat top is designated as the downward region if both the left and right sides of the flat top are within the downward region;

• and the remaining flat tops, identified as the complex region, have their trends disregarded.

Here, we calculate the global trend using equation (13).

$$GlobalTrend=\left\{\begin{array}{cc}{\displaystyle Upward,}\hfill & \mathrm{if}\frac{\mathrm{Length}\left(UpwardRegion\right)}{\mathrm{Length}(UpwardRegion+DownwardRegion)}>0.60\hfill \\ {\displaystyle Downward,}\hfill & \mathrm{if}\frac{\mathrm{Length}\left(DownwardRegion\right)}{\mathrm{Length}(UpwardRegion+DownwardRegion)}>0.60\hfill \\ {\displaystyle Unclear,}\hfill & \mathrm{Others}\hfill \end{array}\right.$$

(13)

In equation (13), the global trend is represented as a ratio, which is subsequently assessed by establishing a threshold value. The selection of this threshold bears a direct influence on the efficacy of the global trend determination. We used the 22 days of lightning weather data analysed statistically in the next section and set the threshold to 0.6 by the 95th percentile method. The outcomes of the trend computations are illustrated in Figure 3. The corrected upward (downward) trend range is depicted in green (red). By employing equation (13), the signal in Figure 3 is characterized as a global uptrend. This methodology provides a more holistic examination of the signal trend, as opposed to solely relying on slope for the determination of the global trend. Furthermore, the presence of a data window of size s (equation (8)) imparts partial resistance to interference, ensuring that substantial local signal fluctuations do not compromise the global trend determination.

Figure 3. Global trend of the Morpho with rectification.

Figure 3. Global trend of the Morpho with rectification.

3. Time-frequency spectrum feature statistics

Figure 4. The AEF changes during a lightning event.

Figure 4. The AEF changes during a lightning event.

Figure 5. Time-frequency spectrum of AEF by WT for a lightning event.

Figure 5. Time-frequency spectrum of AEF by WT for a lightning event.

Table 1. Spectral Bandwidth and STD statistical results of lightning and non-lightning processes.

Table 2. Statistical results of energy differences between warning and de-warning of the original AEF data, CEEMDAN-SG denoised data, and EEWT-ASG denoised data.

Figure 6. The AEF energy changes during a lightning process.

Figure 6. The AEF energy changes during a lightning process.

4. Lightning risk warning method and evaluation

4.1. Lightning risk warning method based on AEF signal

The lightning risk warning method proposed in this paper consists of two components: warning and de-warning. The "global trend" referenced here is computed using the method detailed in Section 2.2. The criteria for warning and de-warning are outlined below.

Warning conditions:

• $B\left(n\right)$ and $E^{{}^{\prime }}\left(n\right)$ of the AEF signal show a global upward trend over 20 min;

• $|\Delta B_{\mathrm{avg}}|$> 0.5, $E_{\mathrm{avg}}^{{}^{\prime }}$> -20 dB, and $STD$> 0.5 in 10 min;

During the design of de-warning conditions, we found that relying solely on AEF energy at the moments of warning and de-warning led to inaccuracies in identifying the conclusion of a minor fraction of lightning events. To guarantee the smooth operation of the entire method, we devised an alternative Scheme II as a supplementary measure for the lightning risk warning method in cases where the de-warning signal is unable to be issued properly.

De-warning conditions:

Scheme I

• $B\left(n\right)$ and $E^{{}^{\prime }}\left(n\right)$ of the AEF signal show a global downward trend over 20 min;

• $|E_{\mathrm{avg}}^{{}^{\prime }}-E_{\mathrm{w}}^{{}^{\prime }}|$≤ 1dB and $STD$≤ 1 in 10 min;

Scheme II

• The sum of the 10 judgments for $|\Delta B_{\mathrm{avg}}|$ was less than 0.1;

• $STD$≤ 1 in 10 min;

In the Scheme I, $E_{\mathrm{w}}^{{}^{\prime }}$ is the minimum value of $E_{\mathrm{avg}}^{{}^{\prime }}$ among the 10 judgments before the lightning risk warning method is judged as a warning. $E_{\mathrm{w}}^{{}^{\prime }}$ changed with the value of $E_{\mathrm{avg}}^{{}^{\prime }}$ each time the lightning risk warning method was judged to be a warning state. A flowchart of the lightning warning method is presented in Figure 7.

Figure 7. Flow chart of the lightning risk warning method.

Figure 7. Flow chart of the lightning risk warning method.

4.2. Evaluation metrics

In accordance with the scholarly definitions of the Area of Concern (AOC) and Warning Area (WA) within the "Two Area Method" [7,26,27], we establish the AOC and WA as illustrated in Figure 8. Here, the circle’s center represents the location of the AEFM, while the AOC encompasses the area surrounding the AEFM that requires protection from lightning hazards. Additionally, the WA constitutes the outer region encircling the AOC.

To assess the performance of the warning behavior, we utilized five metrics: Probability of Detection (POD), Miss Alarm Rate (MAR), False Alarm Rate (FAR), Critical Success Index (CSI), and Mean Warning Lead Time (WLT). These metrics are represented by equations (16), (17), (18), (19), and (20), respectively.

$$POD=\frac{EA}{EA+FTW}$$

(16)

$$MAR=\frac{FTW}{EA+FTW}=1-POD$$

(17)

$$FAR=\frac{FA}{FA+EA}$$

(18)

$$CSI=\frac{EA}{FA+EA+FTW}$$

(19)

$$WLT=\frac{{\sum }_{i=1}^{EA}(tim{e}_i^{{}^{\prime }}-tim{e}_i)}{EA}$$

(20)

In these metrics, $EA$ represents the number of effective alarms, referring to warnings issued prior to the first cloud-to-ground (CG) flash within a warning cycle. Conversely, $FTW$ denotes the number of failures to warn, which occur when warnings are either issued or absent after the first CG flash within a warning cycle. Additionally, $FA$ signifies the number of false alarms, indicating the absence of a CG flash during a warning cycle. Lastly, $tim{e}_i^{{}^{\prime }}$ and $tim{e}_i$ correspond to the time of the first CG flash as detected by the Guangdong-Hong Kong-Macau Lightning Location System (GHMLLS), and the warning time provided by the AEF-based lightning risk warning method, respectively.

Figure 8. Configuration of AOC and WA.

Figure 8. Configuration of AOC and WA.

4.3. Results and Analysis

We evaluated the lightning risk warning method proposed in this paper using AEF data with 83 lightning events, employing the evaluation metrics in Section 4.2. In addition, using 123 days of AEF data for non-lightning events, we employed POD to evaluate the performance of the lightning risk warning method for non-lightning events. The statistical results are presented in Table 3.

Table 3. Statistics on the performance of different lightning risk warning methods.

Table 3. Statistics on the performance of different lightning risk warning methods.

Method	Lightning					Non-lightning
	POD	MAR	FAR	CSI	WLT/min	POD
EFAI [28]	18.07%	81.93%	80.52%	0.11	19.19	42.61%
EFDI [28]	91.57%	8.43%	60.00%	0.38	32.48	80.00%
Lu [13]	18.07%	81.93%	70.59%	0.13	10.04	98.26%
We proposed	77.11%	22.89%	40.19%	0.51	22.27	90.24%

In Table 3, we compare our proposed method with previous studies in the literature [13,28]. It was found that the EFDI method had the highest Probability of Detection (POD) during lightning events, but its False Alarm Rate (FAR) reached 60.00%, indicating low reliability [29]. Lu’s method had the highest POD of 98.26% during non-lightning events, but it performed poorly with a POD of 18.07% during lightning events and failed to meet the warning goal. Our proposed method achieved the lowest FAR (40.19%) and highest Critical Success Index (CSI) (0.51) among the four methods during lightning events, while maintaining a reasonable POD. Additionally, our method also had a high POD (90.24%) during non-lightning events, although it was lower than that of Lu’s method. We also visualized the warning lead time of the warning methods in this paper (Figure 9), and it can be seen that the warning lead times are mostly within the range of 0 to 20 minutes.

Figure 9. The distribution of warning lead time.

Figure 9. The distribution of warning lead time.

5. Discussions and conclusions

References

1. Holle, R.L. A summary of recent national-scale lightning fatality studies. Weather, Climate, and Society 2016, 8, 35–42. [Google Scholar] [CrossRef]
2. Yin, Q.; Liu, H.; Fan, X.; Zhang, Y.; Zhuang, Y.; Wang, F.; Du, H.; Huang, X.; Chen, S.; Chen, L. Lightning fatalities in China, 2009-2018. Journal of Agricultural Meteorology 2021, 77, 150–159. [Google Scholar]
3. Bernard, M.; Underwood, S.J.; Berti, M.; Simoni, A.; Gregoretti, C. Observations of the atmospheric electric field preceding intense rainfall events in the Dolomite Alps near Cortina d’Ampezzo, Italy. Meteorology and Atmospheric Physics 2020, 132, 99–111. [Google Scholar] [CrossRef]
4. Odzimek, A.; Baranski, P.; Kubicki, M.; Jasinkiewicz, D. Electrical signatures of Nimbostratus and Stratus clouds in ground-level vertical atmospheric electric field and current density at mid-latitude station Swider, Poland. Atmospheric Research 2018, 209, 188–203. [Google Scholar] [CrossRef]
5. Hoeft, R.; Wakefield, C. Evaluation of the electric field mill as an effective and efficient means of lightning detection. International Aerospace and Ground Conference on Lightning and Static Electricity, National Interagency Coordinating Group, Atlantic City, NJ, US, 1992.
6. Montanya, J.; Bergas, J.; Hermoso, B. Electric field measurements at ground level as a basis for lightning hazard warning. Journal of Electrostatics 2004, 60, 241–246. [Google Scholar] [CrossRef]
7. Murphy, M.J.; Holle, R.L.; Demetriades, N.W. Cloud-to-ground lightning warnings using electric field mill and lightning observations. 20th international lightning detection conference, 2008, pp. 21–23.
8. Aranguren, D.; Montanya, J.; Sola, G.; March, V.; Romero, D.; Torres, H. On the lightning hazard warning using electrostatic field: Analysis of summer thunderstorms in Spain. Journal of Electrostatics 2009, 67, 507–512. [Google Scholar] [CrossRef]
9. Srivastava, A.; Mishra, M.; Kumar, M. Lightning alarm system using stochastic modelling. Natural Hazards 2015, 75, 1–11. [Google Scholar] [CrossRef]
10. Ying, L.I.; Wang, Z.; XIAO. FFT analysis on atmospheric electric field and the applications in thunderstorm prediction. Journal of the Meteorological Sciences 2013, 33, 66–70. [Google Scholar]
11. Xu, D.; Wang, Z.; Zeng, Q.; Xue, A. EMD Based Analysis of Atmospheric Electric Field Data. Meteorological Science and Technology 2013. [Google Scholar]
12. Kang, H.; Liu, C.; Jiang, X. Weather Recognition Algorithm Based on the Characteristics of Atmospheric Electric Field Signal [J]. Computer Simulation 2014, 31, 312–315. [Google Scholar]
13. Lu, Y.; Zhou, Z.; Gu, S.; Wu, D.; Guo, J.; Tao, H. Research on lightning warning method based on the characteristics of atmospheric electric field. 2016 33rd International Conference on Lightning Protection (ICLP). IEEE, 2016, pp. 1–5.
14. Yang, X.; Xing, H.; Zhuang, L. A thunderstorm cloud point charge localization method based on CEEMDAN and SG filtering. IEEE Access 2021, 9, 17049–17059. [Google Scholar] [CrossRef]
15. Torres, M.E.; Colominas, M.A.; Schlotthauer, G.; Flandrin, P. A complete ensemble empirical mode decomposition with adaptive noise. 2011 IEEE international conference on acoustics, speech and signal processing (ICASSP). IEEE, 2011, pp. 4144–4147.
16. Gilles, J. Empirical wavelet transform. IEEE transactions on signal processing 2013, 61, 3999–4010. [Google Scholar] [CrossRef]
17. Gilles, J.; Heal, K. A parameterless scale-space approach to find meaningful modes in histograms—Application to image and spectrum segmentation. International Journal of Wavelets, Multiresolution and Information Processing 2014, 12, 1450044. [Google Scholar] [CrossRef]
18. Hu, Y.; Li, F.; Li, H.; Liu, C. An enhanced empirical wavelet transform for noisy and non-stationary signal processing. Digital signal processing 2017, 60, 220–229. [Google Scholar] [CrossRef]
19. Lee Rodgers, J.; Nicewander, W.A. Thirteen ways to look at the correlation coefficient. The American Statistician 1988, 42, 59–66. [Google Scholar] [CrossRef]
20. John, A.; Sadasivan, J.; Seelamantula, C.S. Adaptive Savitzky-Golay filtering in non-Gaussian noise. IEEE Transactions on Signal Processing 2021, 69, 5021–5036. [Google Scholar] [CrossRef]
21. Muthusivagami, R. A survey on ECG signal denoising using S-transform and SG filtering. International Journal of Advanced Research 2020, 8, 332–336. [Google Scholar] [CrossRef]
22. Gilles, J.; Tran, G.; Osher, S. 2D empirical transforms. Wavelets, ridgelets, and curvelets revisited. SIAM Journal on Imaging Sciences 2014, 7, 157–186. [Google Scholar] [CrossRef]
23. Yang, P.; Peng, C.; Fang, D.; Wen, X.; Xia, S. Design, fabrication and application of an SOI-based resonant electric field microsensor with coplanar comb-shaped electrodes. Journal of Micromechanics and Microengineering 2013, 23, 055002. [Google Scholar] [CrossRef]
24. Yang, P.; Chen, B.; Wen, X.; Peng, C.; Xia, S.; Hao, Y. A novel MEMS chip-based atmospheric electric field sensor for lightning hazard warning applications. 2015 IEEE SENSORS. IEEE, 2015, pp. 1–4.
25. Zhang, Y.; Lu, W.; Chen, L.; Wu, B.; Qi, Q.; Ma, Y.; Zhang, Y.; Zheng, D.; Yan, X.; Meng, Q. Evaluation of GHMLLS Performance Characteristics Based on Observations of Artificially Triggered Lightning. Journal of Applied Meteorological Science 2022, 33, 329–340. [Google Scholar]
26. Zhou, J.; Zeng, Q.; Shah, A.B.; Jiao, X. A lighting warning algorithm using electric filed mills network and lightning locating system. 2016 12th International Conference on Natural Computation, Fuzzy Systems and Knowledge Discovery (ICNC-FSKD). IEEE, 2016, pp. 1631–1635.
27. Li, X.; Hayawi, K.; Chen, Y.; Chang, S.Y.; Wen, H.; Ho, P.H.; Yang, L.; Yin, Q. False Data Injection Attack on Atmospheric Electric Field in Thunderstorm Warning. 2022 International Conference on Computing, Communication, Perception and Quantum Technology (CCPQT). IEEE, 2022, pp. 219–223.
28. Zeng, Q.; Wang, Z.; Guo, F.; Feng, M.; Zhou, S.; Wang, H.; Xu, D. The application of lightning forecasting based on surface electrostatic field observations and radar data. Journal of Electrostatics 2013, 71, 6–13. [Google Scholar] [CrossRef]
29. Cooper, M.A.; Holle, R.L. Reducing lightning injuries worldwide; Springer, 2019.