Automatic Recognition of Constant-Frequency Electromagnetic Disturbances Observed by the Electric Field Detector on Board the CSES

Abstract

Since the CSES (China Seismo-Electromagnetic Satellite) has been in orbit, it has detected a large number of constant-frequency electromagnetic disturbances (CFEDs), which are horizontal lines on the spectrum. In this paper, we present an algorithm for automatic recognition of CFEDs based on computer vision technology. The relevant results are of great significance for analysis of perturbation events and mining of the transformation laws of global space events. First, a grayscale spectrogram is obtained; then, a horizontal convolution kernel is used to enhance the horizontal edge features of the grayscale graph, and finally, black-and-white binarization is performed to complete data preprocessing. The preprocessed data are then fed into an unsupervised cluster model for training and recognition to realize automatic recognition of CFEDs. Experimental results show that the CFED recognition algorithm proposed in this paper is effective, with a recognition accuracy of more than 98%.

1. Introduction

Since the 1980s, electromagnetic satellites have captured a large number of electromagnetic disturbance phenomena [1,2,3,4,5,6,7] caused by solar magnetic storms, substorms, lightning, atmospheric conditions, tides, artificial VLF transmitter stations, power systems, satellite platforms, etc. Different types of phenomena appear in different shapes in their spectrogram [8]. For example, electromagnetic disturbances produced by the artificial VLF transmitting station, PHRL generated by power systems, and disturbances of the satellite platform itself could cause changes in the spatial physics environment, all of which present as a horizontal shape on the spectrogram. The artificial VLF transmitting station transmits continuous electromagnetic waves at a constant frequency in the range of 10 kHz to 50 kHz, causing a horizontal shape on the spectrogram [9]. Furthermore, the PHRL recorded by the satellite can display a series of parallel lines in the electromagnetic spectrogram. These lines have constant frequency intervals of 50/100 Hz or 60/120 Hz. PHRL is typically distributed in the extremely low frequency (ELF) or VLF range [8,10,11].

In February 2018, the CSES was successfully launched an has since been in orbit for more than 4 years and generated a vast body of data. Faced with the challenge of a massive amount of data, it is urgent to automatically recognize abnormal constant-frequency electromagnetic disturbances (CFEDs). Recently, some progress has been made in automatic recognition of CFEDs. This work can be divided into two categories. The first group is based on the recognition of CFEDs by Hough transformation [12]. This method detects line segments in many directions, although it is time-consuming and leads to unstable performance. The second group is based on density statistics [13], employing a horizontal convolution kernel to enhance the features of the horizontal edge, giving rise to improved recognition performance. However, the robustness of these methods is relatively poor. To overcome this limitation, we performed cluster analysis to propose a novel CFED recognition method.

Cluster analysis is a kind of unsupervised machine learning algorithm that groups objects such that those in the same cluster are similar to each other [14]. Compared with deep learning methods, clustering algorithms not only do not require labels and datasets for training but are also more explanatory [15]. Recently, cluster analysis algorithms have been widely used in geophysics and other related fields. For example, the authors of [16] used cluster analysis to classify types of electrical sounding curves. Among many types of clustering algorithms, the most popular is K-means, which is easy to implement, with a moderate running time [17,18,19,20]. Based on K-means, we propose a novel CFED recognition algorithm to be applied to spectrograms from EFD data. The proposed algorithm provides innovations: (1) an error rate of 0 and (2) robust performance.

2. Data Collection

The CSES is China’s first geophysical exploration satellite. It is mainly used to acquire global electromagnetic field, plasma, high-energy particle observation and other data, are used to provide services such as short-term seismic prediction and geophysical research [21,22,23,24]. The CSES has an orbital inclination of 97.4°, which is a solar-synchronous orbit. The regression period is 5 days. It takes approximately 94 min to complete an orbit around the Earth. Observation data are stored as ascending track (night) and descending track (day), and an observation of each half track (ascending track/descending track) is approximately 34 min. The spatial resolution of adjacent ascending (or descending) tracks within the same day is approximately 2000 km. The CSES carries eight scientific payloads [25,26,27,28,29]: an induction magnetometer, a high-precision magnetometer [30], an electric field detector (EFD), a GNSS occultation receiver [31], a plasma analyzer [32,33], a high-energy particle detector [34], a Langmuir probe [35] and a tri-frequency beacon transmitter [22]. The electric field detector (EFD) detects electric fields in space. Its detection frequency band is divided into ULF (0~16 Hz), ELF (6 Hz~2.2 kHz) and VLF (1.8 kHz~35 MHz). It can not only provide basic data for the study of Sun–Earth space physics; space weather; and interactions between the ionosphere and the upper atmosphere, magnetosphere and other related spheres and their effects but also provides data application services for seismic observation and research [36].

The EFD waveform data of the VLF band-Z component in this referenced in this paper were obtained by the CSES. A spectrogram is obtained by short-time Fourier transformation (STFT) of waveform data. In Figure 1, the y-axis represents the frequency range; the x-axis represents time, longitude and latitude of the satellite orbit; and the color bar on the right represents the intensity of electromagnetic waves.

3. CFED Automatic Recognition Algorithm

The goal of the algorithm is to recognize horizontal lines on the spectrogram; as shown in Figure 1, several horizontal lines are observed easily. The framework of the proposed algorithm is provided in Figure 2. The solid line represents the training flow, and the dotted line is the test flow. By following the training flow, a recognition model for CFEDs is obtained, and by following the test flow, the recognition results are obtained.

3.1. Grayscale Processing

Information about horizontal shapes in CFEDs is more important than color information. Therefore, we first process the time-frequency diagram from RGB (red, green and blue) to gray. The gray image is obtained using the red, blue or green channel in RGB color space, i.e., ImgG = RGB.R [13]. Note that when the color bar of the spectrogram is changed, the grayscale method is changed accordingly. Experiments show that when the color bar ranges from −10 to 12, the best result of the graying experiment is achieved using the red channel. When the color bar range is −2 to 1, the best result is achieved using the blue channel.

3.2. Horizontal Feature Enhancement

To enhance CFED features for easy recognition, a convolution kernel [13] is used to process the ImgG image. The kernel is provided by Equation (1).

(1)$\mathrm{kernel}=[1,0,-1]$

The method is expressed as Equation (2).

(2)$\mathrm{dst}(x,y)={\displaystyle \sum _{\begin{array}{l}0\le x^{\prime }<\mathrm{kernel}.\mathrm{cols},\\ 0\le y^{\prime }<\mathrm{kernel}.\mathrm{rows}\end{array}}\mathrm{kernel}(x^{\prime },y^{\prime })\times \mathrm{src}(x+x^{\prime }-\mathrm{anchor}.x,y+y^{\prime }-\mathrm{archor}.y)}$

where src is the original image; dst is the target image; kernel is the convolution kernel; anchor is the anchor point of the kernel, indicating the relative position of the kernel filter point; x and y are the coordinates of pixels on the original graph; and x′ and y′ are the coordinates of points in the convolution kernel.

Figure 3a shows the effect of a 5 × 5 pixel spectrogram enlarged 500 times without distortion. A horizontal line is clearly visible above the background. The operational process is shown in Figure 3b, where src is the red channel value of Figure 3a, and dst is the result of the convolution operation. The results show that the convolution operation greatly enhances the edge features of the horizontal line.

3.3. Binarization

In order to separate the object from the background more quickly and accurately, black-and-white binarization is carried out on the convolution result so that the whole image presents an obvious black-and-white effect. Binarization is expressed by Equation (3).

(3)$\mathrm{dstmap}(px,py)=\{\begin{array}{ll}maxVal& if\mathrm{srcmap}(px,py)>tresh\\ 0& \mathrm{otherwise}\end{array}$

where srcmap is the target image after convolution; px and py are the coordinates of a pixel on the target image; maxVal represents the maximum number of pixels, which is set to 255 in this paper; and thresh is the threshold value, which is set to 10 in this paper. if the pixel value of (px,py) is less than thresh, it is 0; otherwise, it is 255. dstmap is the output image after binarization. The results after binarization of Figure 3b (dst) are shown in Figure 4. After binarization, only the pixel value on the line is 255, and the rest is 0.

3.4. K-Means Clustering

The preprocessed data described above are fed into the K-means model. The steps are as follows.

Step 1: Randomly select k sample data points from the sample data to be clustered as the initial clustering center {µ1,µ2,…,µn} (k = 2).

Step 2: In the K-means model, if the sum of the distance between sample points in the cluster and the center of the cluster is small, the samples in the cluster are similar, and the clustering effect is better. The calculation formula is expressed as Equation (4):

(4)$d={\displaystyle \sum _{i=1}^k{\displaystyle \sum _{p\in C_i}{\Vert p-{\mu }_i\Vert }^2}}$

where ${\mu }_{\mathrm{i}}$ is the central sample of cluster $C_{\mathrm{i}}$, and $p$ is the sample point.

Step 3: Recalculate the clustering centers for the newly obtained k clusters. The formula for calculating the center of the cluster is expressed as Equation (5):

(5)${\mu }_{\mathrm{i}}=\frac{1}{C_{\mathrm{i}}}{\displaystyle \sum _{p\in C_i}p}$

Repeat steps 2 and 3 until the termination condition ($\left|{\mu }_{\mathrm{n}+1}-{\mu }_{\mathrm{n}}\right|\le \epsilon$) is satisfied or the preset number of iterations is reached; then, the algorithm terminates.

In this paper, the number of iterations is set to 10, and the iteration with the lowest d value is considered the clustering result.

4. Experiment and Analysis

4.1. Experimental Methods

Python 3.7 is used to implement the detection process described above using the CV2 image library and the SKLearn clustering algorithm library. The following experimental scheme was formulated.

• (1). Data collection

In this experiment, we collected more than 500 spectrograms from EFD data recorded in April 2020, each of which contains several horizontal lines.

• (2). Data preprocessing

First, the spectrogram was processed by taking data from the red channel of the RGB image to generate the gray image; then, the horizontal edge features were extracted, and finally, they were fed into the model after binarization. Some of the results obtained by the above methods are shown in Figure 5.

• (3). Model training

There are three commonly used unsupervised clustering models: BIRCH, K-means, and DBSCAN.

We selected spectrograms with clean CFEDs and processed them to generate binarized data. The binarized data were then fed into three training models.

• (4). Horizontal line recognition

The preprocessed datasets were fed into the three models for line identification, and the experimental results were counted and analyzed to determine which type of model to choose.

4.2. Optimal Clustering Coefficient Selection

Taking the spectrograms in Figure 5 as an example, binarization data can be visualized according to the line graph to check the distribution characteristics of the data. As shown in Figure 6, the x-axis represents the number of rows, and the y-axis represents the sum of pixel values for each row. In the original spectrogram shown in Figure 5a, there are five clearly visible horizontal lines and five sum values greater than 30,000, corresponding to the data distribution characteristics shown in Figure 6.

The DBSCAN (density-based spatial clustering) algorithm uses spatial index technology to search the neighborhood of an object. It treats a cluster as a dense object area separated by a low-density area in the data space. The binarized dataset does not fit the DBSCAN model usage scenario. We also experimentally verified that the DBSCAN classification effect is very poor, and it can hardly recognize horizontal lines. Two other BIRCH and K-means models are discussed below.

BIRCH (balanced iterative reduction and clustering using hierarchies) is a comprehensive hierarchical clustering algorithm. It uses the concepts of clustering feature (CF) and clustering feature tree (CF Tree) to summarize the description of clustering.

The clustering coefficient (k) is an important parameter of the clustering algorithm. With the exception of k, all other parameters in the algorithm adopt default values. There are two common ways to select k values.

• (1). Practical application

The goal is to recognize horizontal lines. According to this practical application requirement, the clustering coefficient (k) should be set to 2; that is, the target is divided into two categories: pixel points on the line and pixel points outside the line.

• (2). Silhouette method

The silhouette method considers the density and dispersion of clusters comprehensively. If the data are divided into ideal k clusters, the samples within the corresponding cluster are very dense, whereas the samples between clusters are very scattered. The silhouette value (s(i)) is expressed by Equation (6). For the ith sample, the larger the s(i) value, the more reasonable the sample distribution is.

(6)$s(i)=\frac{b(i)-a(i)}{\max \{b(i),a(i)\}}$

where a(i) represents the average distance between sample i and other sample points in the same cluster. A low a(i) value indicates that sample i should be clustered in this cluster. The mean distance (bij) between sample i and all samples of other clusters (Cj) is called the dissimilarity between sample i and cluster Cj, indicating the dissimilarity between clusters of sample i: b(i) = min{bi1,bi2… bik}. The silhouette coefficient is the mean of the silhouette coefficient of all sample points. Figure 7 shows the change in silhouette coefficient when the k value changes. When k = 2, the silhouette coefficient value is the maximum, indicating that k = 2 is the best optional reference value.

4.3. Experimental Results

After the clustering operation, the data sample is divided into two clusters, labeled 0 or 1, namely linear cluster and nonlinear clusters, respectively. Then, each row on the original spectrogram is traversed line by a line. If it belongs to the row cluster, the pixel points obtained after binarization are marked in red. Experimental results from both models are shown in Figure 8. According to the recognition results, K-means is superior to BIRCH.

(a)

(b)

In the statistical analysis of detection results, if some horizontal lines are intermittent or indistinct, they are not included in the statistical scope of the experiment. However, if the experiment is recognized by the algorithm, it is not a detection error; if it is not recognized, it is not a missed detection. Statistics of experimental results are shown in Table 1.

Therefore, the K-means algorithm with k = 2 is selected for actual line recognition.

4.4. Comparison with Existing Algorithms

At home and abroad, the recognition technology of CFED in the spectrogram is mainly based on Hough line transform [12]. This recognition technology is aimed at straight lines in various directions, and the recognition efficiency and accuracy are relatively low. The authors of [13] proposed an improved row density statistics algorithm, which greatly improved the efficiency and accuracy of recognition. However, the robustness of this method is poor, and the threshold must be set manually. If the threshold is too high, the default rate increases, and if the threshold is too low, many lines are incorrectly recognized as straight lines. The K-means algorithm based on horizontal a convolution kernel proposed in this paper not only has a high accuracy rate and almost zero error rate achieves robust performance. The experimental statistical results of the three methods are shown in Table 2.

5. Discussion

In this section, the influence of horizontal line recognition results is analyzed and discussed in depth, and we explore future research prospects for horizontal line recognition.

5.1. Cause Analysis of Missing and Incorrect Detection

The reasons for missing or misidentified horizontal straight lines are discussed below.

First, the CFED signal itself is relatively weak in the faint visible spectrum, as shown in Figure 9a; second, the horizontal line is disturbed by other spatial electromagnetic waves, and its shape on the spectrogram changes as a result of contrast reduction or clarity reduction, as shown in Figure 9b. Third, the CFED signal itself is interrupted, as shown in Figure 9c.

5.2. Extended Discussion

According to electromagnetic wave propagation theory, most of the energy in VLF/LF wave propagation is confined to the waveguide layer. When the low ionosphere over the propagation path is disturbed by various factors, such as solar flare, magnetic storms, lightning flash discharge, earthquakes, etc., VLF/LF artificial source signal anomalies occur [37]. According to the linear detection algorithm proposed in this paper, when some CFEDs are disturbed by other spatial electromagnetic waves, the CFED signal is abnormal, resulting in the weakening, interruption, or even disappearance of the horizontal straight line on the spectrum diagram, and the corresponding detection results can also be reflected by the marks. The interruption of the line or the weakening of the signal is helpful for the detection of spatial electromagnetic anomaly disturbances. The next step is to analyze the influence of electromagnetic disturbances caused by the 2020 Cuban earthquake on CFEDs.

At 19:10 on 28 January 2020, an earthquake occurred 125 km northwest of the port city of Lusi, Jamaica, at a depth of 10 km (Figure 10). The ionospheric disturbances over the seismogenic zone in the southern seas of Cuba before the earthquake were analyzed based on the observational data of the VLF frequency band of the space electric field.

The geomagnetic activity index is closely related to ionospheric changes in space. In order to obtain and identify anomalous information of ionospheric disturbances related to earthquakes, it is necessary to exclude disturbances caused by known geomagnetic activities. The ∑Kp and Dst indices in the two months prior to the earthquake were collected and plotted; the changes in activity are shown in Figure 11, with the study period of this earthquake covering the geomagnetic quiet period [26,38].

By using the automatic CFED detection and labeling algorithm proposed in this paper, abnormal electromagnetic disturbance in the earthquake area can be judged by statistical point density and sparse morphological changes, as shown in Figure 12.

6. Conclusions

In this paper, we propose a method for automatic recognition of CFEDs on the spectrogram. Three types of clustering models were trained and fed into the data prediction model. Through experimental statistical analysis and various evaluation methods, we found that K-means with k = 2 had the best recognition effect on horizontal lines, with an accuracy rate of more than 98%. It was helpful to further excavate abnormal disturbance events and study temporal and spatial variation of the space ionization field.

However, because some horizontal electromagnetic signals are affected by other spatial electromagnetic interference factors, horizontal lines become weak or disappear, resulting in missed detection. On the other hand, analysis shows that if the mark shape of the detected horizontal line points changes, such as a certain section in the middle becoming sparse, this area can be marked as abnormal, these abnormal areas can be mined, and the causes of the abnormalities can be discussed and analyzed.

Accordingly, future work will focus mainly on the following two areas: (1) excavation and analysis of the phenomena of weakening or disappearance of CFED caused by other electromagnetic interference and (2) automatic recognition of other types of spatial electromagnetic interference events using computer vision processing based on the scheme proposed herein.

Author Contributions

Conceptualization and methodology, Y.H.; algorithm implementation, J.Y.; data analysis and conclusion, Q.O., J.H., Z.L.; software and investigation, Y.Z., Y.W.; writing—review and editing, X.S. and Z.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Publicly available datasets were analyzed in this study. The CSES electric field data can be found here: (www.leos.ac.cn, accessed on 1 September 2022).

Acknowledgments

This work made use of the data from the CSES mission, a project funded by the China National Space Administration (CNSA) and the China Earthquake Administration (CEA). Thanks to the CSES team for the data (www.leos.ac.cn, accessed on 1 September 2022).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. VLF from 6:40 a.m. to 7:14 a.m. on 25 January 2020 CSES waveform spectrogram.

Figure 2. Recognition algorithm (solid line is training flow; dotted line is testing flow).

Figure 3. Convolution operation: (a) spectrogram; (b) convolution operation.

Figure 4. Binarization.

Figure 5. Data preprocessing: (a) spectrogram; (b) grayscale; (c) extraction of horizontal features; (d) binarization.

Figure 6. Distribution of data characteristics.

Figure 7. Contour coefficients corresponding to different clustering coefficients in the two clustering algorithms.

Figure 8. Marking of prediction results of different clustering algorithms (k = 2): (a) K-means; (b) BIRCH.

Figure 9. Marking and analysis of test results. (a) weak signal; (b) other electromagnetic wave interference; (c) signal interruption.

Figure 10. Satellite trajectory projection on 28 January 2020.

Figure 11. Changes in the ∑Kp and Dst indices during the selected study period of the sample earthquake.

Figure 12. Automatic detection and labeling results.

Table 1

Experimental results.

	Accuracy Rate	Missed Rate	Error Rate
K-means	0.989 ± 0.02	0.011 ± 0.03	0
BIRCH	0.8712 ± 0.13	0.1388 ± 0.06	0

Table 2

Experimental results.

	Accuracy Rate	Miss Rate	Error Rate
Hough	91.12 ± 0.12	8.88 ± 0.18	15.38 ± 0.01
Row density statistics	99.12 ± 0.01	0.0088 ± 0.01	8.1 ± 0.03
K-means	0.989 ± 0.02	0.011 ± 0.03	0

References

1. Yan, R.; Shen, X.; Huang, J.; Wang, Q.; Chu, W.; Liu, D.; Yang, Y.; Lu, H.; Xu, S. Examples of unusual ionospheric observations by the CSES prior to earthquakes. Earth Planet. Phys.; 2018; 2, pp. 79-90. [DOI: https://dx.doi.org/10.26464/epp2018050]

2. Larkina, V.I.; Nalivaiko, A.V.; Gershenzon, N.I.; Gokhberg, M.B.; Liperovskii, V.A.; Shalimov, S.L. Intercosmos-19 observations of VLF emissions associated with seismic activity. Geomagn. Aeron.; 1983; 23, pp. 842-846.

3. Chmyrev, V.M.; Isaev, N.V.; Bilichenko, S.V.; Stanev, G. Observation by space-borne detectors of electric fields and hydromagnetic waves in the ionosphere over an earthquake center. Phys. Earth Planet. Inter.; 1989; 57, pp. 110-114. [DOI: https://dx.doi.org/10.1016/0031-9201(89)90220-3]

4. Parrot, M.; Mogilevsky, M. VLF emissions associated with earthquakes and observed in the ionosphere and magnetosphere. Phys. Earth Planet. Inter.; 1989; 57, pp. 86-99. [DOI: https://dx.doi.org/10.1016/0031-9201(89)90218-5]

5. Pulinets, S.A.; Legen’ka, A.D. Spatial–Temporal Characteristics of Large Scale Disturbances of Electron Density Observed in the Ionospheric F-Region before Strong earthquakes. Cosm. Res.; 2003; 41, pp. 221-230. [DOI: https://dx.doi.org/10.1023/A:1024046814173]

6. Liu, Y.M.; Wang, J.S.; Xiao, Z.; Suo, Y. A Possible Mechanism of Typhoon Effects on the Ionospheric F2 layer. Chin. J. Space Sci.; 2006; 26, pp. 92-97.

7. Cai, J.T.; Zhao, G.Z.; Zhan, Y.; Tang, J.; Chen, X.B. The study on ionospheric disturbances during earthquakes. Prog. Geophys.; 2007; 22, pp. 695-701.

8. Parrot, M. DEMETE Robservations of manmade waves that propagate in the ionosphere. Comptes Rendus Phys.; 2018; 19, pp. 26-35. [DOI: https://dx.doi.org/10.1016/j.crhy.2018.02.001]

9. Ryu, K.; Lee, E.; Chae, J.; Parrot, M.; Pulinets, S. Seismo-ionospheric coupling appearing as equatorial electron density enhancements observed via DEMETER electron density measurements. J. Geophys. Res. Space Phys.; 2015; 119, pp. 8524-8542. [DOI: https://dx.doi.org/10.1002/2014JA020284]

10. Zhima, Z. Electromagnetic Waves Induced by Strong Earthquake and Non-Earthquake Sources in Geospace. Doctoral Dissertation; Beijing University of Aeronautics and Astronautics: Beijing, China, 2014.

11. Zhang, X.M.; Qian, J.D.; Ying, Y.; Zhao, S.; Lou, W. Space Electromagnetic Disturbances Induced by High-voltage Power Lines. Earthquake; 2019; 39, 11.

12. Han, Y.; Yuan, J.; Feng, J.; Yang, D.; Huang, J.; Wang, Q.; Shen, X.; Zeren, Z. Automatic detection of “horizontal” electromagnetic wave disturbance in the data of EFD on ZH-1. Prog. Geophys.; 2021; 36, pp. 2303-2311.

13. Han, Y.; Yuan, J.; Feng, J.L.; Yang, D.; Huang, J.; Wang, Q.; Shen, X.; Zeren, Z. Automatic detection of horizontal electromagnetic wave disturbance in EFD data of Zh-1 based on horizontal convolution kernel. Prog. Geophys.; 2022; 37, pp. 11-18.

14. He, L.; Wu, L.D.; Cai, Y. Survey of Clustering Algorithms in Data Mining. Appl. Res. Comput.; 2007; 24, pp. 10-13.

15. Wang, D.; Yuan, S.; Yuan, H.; Zeng, H.; Wang, S. Intelligent velocity picking based on unsupervised clustering with the adaptive threshold constraint. Chin. J. Geophys.; 2021; 64, pp. 1048-1060. (In Chinese)

16. Song, Y.C.; Meng, H.D.; O’grady, M.J.; O’Hare, G.M. Application of cluster analysising geophysical data interpretation. Comput. Geosci.; 2010; 14, 263271. [DOI: https://dx.doi.org/10.1007/s10596-009-9150-1]

17. Huang, Y.; Wang, X.Q.; Zhou, Z.K.; Zou, Z.L.; Xia, G.P.; Weng, A.H. Application of clustering analysis algorithm in three-dimensional magnetotelluric interpretation. Prog. Geophys.; 2019; 34, pp. 568-572. (In Chinese)

18. Hung, C.H.; Chiou, H.M.; Yang, W.N. K-harmonic means data clustering. Appl. Math. Model.; 2013; 37, pp. 10123-10128. [DOI: https://dx.doi.org/10.1016/j.apm.2013.05.052]

19. Macqueen, J. Some methods for classification and analysis of multivariate observations. Proceedings of the 5th Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Statistics; University of California Press: Berkeley, CA, USA, 1967; pp. 281-297.

20. Ahmad, A.; Hashmi, S. K-Harmonic means type clustering algorithm for mixed datasets. Appl. Soft Comput.; 2016; 48, pp. 39-49. [DOI: https://dx.doi.org/10.1016/j.asoc.2016.06.019]

21. Shen, X.; Zong, Q.G.; Zhang, X. Introduction to special section on the China Seismo-Electromagnetic Satellite and preliminary results. Earth Planet. Phys.; 2018; 2, pp. 3-7. [DOI: https://dx.doi.org/10.26464/epp2018041]

22. Zhou, B.; Yang, Y.; Zhang, Y.; Gou, X.; Cheng, B.; Wang, J.; Li, L. Magnetic field data processing methods of the China SeismoElectromagnetic Satellite. Earth Planet. Phys.; 2018; 2, pp. 455-461. [DOI: https://dx.doi.org/10.26464/epp2018043]

23. Cao, J.B.; Yan, C.X.; Lu, L. Non-seismic Induced Electromagnetic Waves in Near-Earth Space. Earthquakes; 2008; 29, pp. 17-25. (In Chinese)

24. Wang, Z.; Zhou, C.; Zhao, S.; Xu, X.; Liu, M.; Liu, Y.; Liao, L.; Shen, X. Numerical Study of Global ELF Electromagnetic Wave Propagation with Respect to Lithosphere–Atmosphere–Ionosphere Coupling. Remote Sens.; 2021; 13, 4107. [DOI: https://dx.doi.org/10.3390/rs13204107]

25. Wang, L.; Hu, Z.; Shen, X. Data processing methods and procedures of CSES satellite. Natl. Remote Sens. Bull.; 2018; 22, pp. 39-55. [DOI: https://dx.doi.org/10.11834/jrs.20187360]

26. Li, Z.; Li, J.; Huang, J.; Yin, H.; Jia, J. Research on Pre-Seismic Feature Recognition of patialElectric Field Data Recorded by CSES. Atmosphere; 2022; 13, 179. [DOI: https://dx.doi.org/10.3390/atmos13020179]

27. Wang, L.; Shen, X.; Zhang, Y.; Zhang, X.; Hu, Z.; Yan, R.; Yuan, S.; Zhu, X. Developing progress of China Seismo-Electromagnetic Satellite project. Acta Seismol. Sin.; 2016; 38, pp. 376–385+509.

28. Wang, X.; Yang, D.; Zhou, Z.; Cui, J.; Zhou, N.; Shen, X. Features of topside ionospheric background over China and its adjacent areas obtained by the ZH-1 satellite. Chin. J. Geophys.; 2021; 64, pp. 391-409.

29. Shen, X.; Zhang, X.; Cui, J.; Zhou, X.; Jiang, W.; Gong, L.; Li, Y.; Liu, Q. Remote sensing application in earthquake science research and geophysical fields exploration satellite mission in China. Natl. Remote Sens. Bull.; 2018; 22, pp. 5-20.

30. Zhou, B.; Cheng, B. Development and calibration of high-precision magnetometer of the China Seismo-Electromagnetic Satellite. Natl. Remote Sens. Bull.; 2018; 22, pp. 64-73.

31. Diego, P.; Huang, J.; Piersanti, M.; Badoni, D.; Zeren, Z.; Yan, R.; Rebustini, G.; Ammendola, A.; Candidi, M.; Guan, Y.-B. et al. The Electric Field Detector on Board the China Seismo Electromagnetic Satellite-In-Orbit Results and Validation. Instruments; 2020; 5, 1. [DOI: https://dx.doi.org/10.3390/instruments5010001]

32. Lu, H.; Shen, X.; Zhao, S.; Liao, L.; Lin, J.; Huang, J.; Zeren, Z.; Sun, F.; Guo, F. Typical event observation of the tri-band beacon onboard the CSES (ZH-1) satellite. Chin. J. Radio Sci.; 2021; 37, pp. 426-433. [DOI: https://dx.doi.org/10.12265/j.cjors.2021106]

33. Liu, C.; Guan, Y.; Zheng, X.; Zhang, A.; Piero, D.; Sun, Y. The technology of space plasma in-situ measurement on the China Seismo-Electromagnetic Satellite. Sci. China Technol. Sci.; 2019; 62, pp. 829-838. [DOI: https://dx.doi.org/10.1007/s11431-018-9345-8]

34. Shen, X.; Zhang, X.; Yuan, S.; Wang, L.; Cao, J.; Huang, J.; Zhu, X.; Piergiorgio, P.; Dai, J. The state-of-the-art of the China Seismo-Electromagnetic Satellite mission. Sci. China Technol. Sci.; 2018; 61, pp. 634-642. [DOI: https://dx.doi.org/10.1007/s11431-018-9242-0]

35. Rui, Y.; Zhe, H.; Lanwei, W.; Yibing, G.; Chao, L. Preliminary data inversion method of Langmuir probe onboard CSES. Acta Seismol. Sin.; 2017; 39, pp. 239-247.

36. Ma, M.; Lei, J.; Li, C.; Zong, C.; Li, S.X.; Liu, Z.; Cui, Y. Design Optimization of Zhangheng-1 Space Electric Field Detector. J. Vac. Sci. Technol.; 2018; 38, pp. 582-589.

37. Zhang, X. The development in seismic application research of VLF/LF radio wave. Acta Seismol. Sin.; 2021; 43, 18.

38. Xu, W. Geomagnetism; Seismological Press: Beijing, China, 2003; pp. 288-296.