EARTHQUAKE PREDICTION EVALUATION BASED ON VLF DATA USING A NOVEL INTERSECTION-UNION METHOD

Electromagnetic phenomena, especially those in the Very Low Frequency/Low Frequency (VLF/LF) bands are promising for shortterm earthquake prediction. Seismo-ionospheric perturbations cause a variety of changes in different receiver-transmitter VLF/LF signal paths. Therefore, independent and simultaneous observations at different points thus in different VLF/LF signal propagation paths are necessary to better predict the earthquake. Most of the previous research on VLF data have been based on one path or limited number of paths which examined perturbations in the time domain and less attention has been paid to estimate the location of the earthquake. In the present research, using wavelet analysis, the temporal variations of seismo-ionospheric perturbations and the approximate time of earthquake are predicted. Clear disturbances are observed two weeks before the Kumamoto earthquake happened in Japan in 2016. The novelty of this study is to present an approach called Intersection-Union method to predict earthquake location. Based on the geometry of a VLF/LF network, the Intersection-Union method was introduced to estimate the earthquake epicenter. This method is based on the overlay of earthquake occurrence probable areas. With simultaneous use of different propagation paths by the Intersection-Union method, an area with a radius of about 300 km was determined as the probable location of the earthquake epicenter. The accuracy of the proposed method is 300 km compared with 1000 km accuracy of other earthquake location prediction scenarios. 1 *Corresponding author


INTRODUCTION
Earthquake is one of the most important natural hazards that has always caused great human and financial damages. Therefore, it is crucial to try to predict earthquakes to reduce physical and human damages and take the necessary measures for crisis management. A number of evidences have been presented in recent decades about earthquake predictors, and especially electromagnetic phenomena as useful factors for short-term earthquake prediction Surkov and Hayakawa, 2014;Hayakawa, 2015). Numerous research have been undertaken on correlation between VLF radio signal perturbations and seismic activities (Rozhnoi et al., 2004;Shvets et al., 2004 a,b;Maekawa et al., 2006;Biagi et al., 2011). As Japan is located in the seismic belt of the world, an integrated network of transmitters and receivers of VLF/LF waves has been established throughout the country. The network has been active for the past two decades, recording VLF/LF data (Ouzounov et al., 2018). Since the VLF/LF method consists of an integrated measurement, it is possible to correlate many pre-earthquake effects to the ionosphere of any earthquakes happening within the sensitive area of the transmitter-receiver great circle. This is the main advantage of VLF/LF method (Hayakawa, 2016). Seismo-ionospheric perturbations are variable in space and time and cause various changes even for an earthquake in a different transmitterreceiver path (Yamauchi et al., 2007). Hence, numerous observations of seismo-ionospheric events in different parts of the world and various transmitter-receiver paths are necessary to understand the related processes, including the mechanism of lithosphere-atmosphere-ionosphere coupling. Research on the interaction between seismic activities and radio wave disturbances has been conducted for decades. Over the past two decades, significant progress has been achieved in the study of earthquake precursors. Among the various types of precursors, there are currently several electromagnetic phenomena that are statistically correlated with earthquakes of a magnitude greater than M= 6 based on long-term data. An obvious example is an ionospheric perturbation, which not only occurs in the lower ionosphere (D region) (Hayakawa et al., 2010;Hayakawa, 2007Hayakawa, , 2011, but also exists in the F region above the ionosphere (Liu, 2009;Liu et al., 2013). Pioneering research has been conducted by Russian and Japanese scientists (Gokhberg et al.,1989;Hayakawa et al., 1996;Molchanov et al., 1998) and a great number of research on the use of the VLF/LF method to study seismo-ionospheric perturbations have been undertaken. The VLF/LF method is becoming a global technique in short-term earthquake prediction, which explains the increasing number of VLF/LF networks in many countries . The main advantage of an integrated VLF measurement is that the observed signal is sensitive to all earthquakes near the propagation path between the transmitter and receiver which facilitates data collection on VLF anomalies . Therefore, a number of analytical studies have been performed on the relationship between VLF/LF anomalies and earthquakes of magnitudes greater than M=6 (Hayakawa et al., 2010;Rozhnoi et al., 2004;Maekawa et al., 2006;Kasahara et al., 2008;Chakrabarti et al., 2010). Wavelet analysis is becoming a common tool for analyzing local power changes within a time series. By converting a time series to frequency-time space, both the dominant modes of variability and the way they change over time can be calculated (Torrence and Compo, 1998). A number of studies have been undertaken on VLF anomalies using wavelet analysis. Righetti et al. (2012) investigated the relationship between LF radio signals received in the European VLF/LF network from July 2009 to April 2011 with earthquakes greater than M=5 using wavelet analysis. Biagi et al. (2019) investigated the seismic activities occurring within the European VLF/LF network from 2016 to 2017 by wavelet analysis. Most of the previous research have been based on one or some limited transmitter-receiver VLF/LF wave paths . Also, most of the previous studies on VLF data have focused on temporal perturbations and less attention has been paid to predict the location of the earthquake (Hayakawa et al., 2010). Therefore, there is a need to develop a method based on several paths simultaneously and to predict the location of the earthquake, which is considered as the main objective of this research. In addition, due to lack of comprehensive reports on VLF/LF data and methods for earthquake prediction from spatial perspective, despite a number of research on the use of geospatial information systems (GIS) in earthquake prediction with various methods, we believe that the present study is among the first researches in using GIS to assist earthquake prediction through VLF/LF signals and data along with its accuracy assessment.

Case Study and Data
The 2016 Kumamoto Earthquake happened in Japan was a huge seismic activity with a sequence of three relatively strong earthquakes on April 14 and 15, 2016 Universal Time (UT). The epicenter of the shock at Kumamoto occurred on April 15, 2016, with a magnitude of M=7.3 and at a shallow epicenter (approximately 10 km), shown as a red dot (EQ) in Figure (1). On April 14, 2016, two strong earthquakes with magnitudes of M= 6.5 and 6.4 have been occurred almost at the same place very close to the epicenter. These two earthquakes are assumed as pre-earthquakes for the mainshock . In this study, Japan VLF/LF network data has been employed. The network consists of 8 VLF receiver stations established throughout Japan illustrated as black dots in Figure  (1). These receivers are Nakashibetsu (NSB), Suttsu (STU), Akita (AKT), Imizu (IMZ), Katsuura (KTU), Kamakura (KMK), Toyohashi (TYH), and Anan (ANA), respectively . The receiver can record signals with a time accuracy of 50 milliseconds to 60 seconds whose sampling with an accuracy of one second was used. Figure (1) shows the relative position of the earthquake epicenter (red dot (EQ)) and the JJI transmitter (VLF) (blue dot) in Kyushu. The path between each VLF observation station and the JJI transmitter is indicated by a black line. In this study, only night time amplitude data were used to improve data quality. The daily amplitude shows very little change for the analysis and is also strongly influenced by sudden ionospheric perturbations (Hayakawa, 2015). The period of night time amplitude is JST 2 = 19 h -29 h or time 10 -20 in UT. The 2016 Kumamoto Earthquakes occurred close to the transmitter, so VLF data was used at all VLF/LF network stations in Japan. Data from seven transmitter-receivers paths were used simultaneously to estimate the location and time of 2 Japan Standard Time the 2016 Kumamoto earthquake in Japan. There are no observations at the NSB station in Hokkaido due to a receiver malfunction during the week before the earthquakes, so NSB observations have not been employed in the analysis.

Figure 1
Relative positions of earthquake epicenter, transmitter, receivers, and receiver-transmitter paths.
A Fresnel zone is an ellipsoidal region of space between and around a transmitter and a receiver. The Fresnel zone is very important in the discussion of the wave sensitive region at VLF/LF frequencies. For ionospheric anomaly detection by VLF waves, the wave-sensitive region for each propagation path is defined by the fifth Fresnel zone (Maekawa et al., 2006;Kasahara et al., 2008). In this research, the fifth Fresnel zone of the receiver-transmitter paths has been used to estimate the location of the earthquake. The fifth Fresnel zones are shown in Figure (2). With the coordinates of seven receivers and one transmitter accessible in the Japan VLF/LF network, the receiver-transmitter paths and their effective area, the spatial layer includes the area covered by the fifth Fresnel zone of the propagation paths and the 300 km area around the transmitter and receivers of VLF/LF network were built in a GIS environment to estimate the location of the earthquake epicenter.

Figure 2
Fifth Fresnel zone of the receiver-transmitter paths.

Methodology
The research methodology is divided into two parts including estimating the location and the time of the earthquake. In the time estimation methods, we measure the anomalies of VLF waves and discover the anomalies at the time of the earthquake. In estimating the earthquake location with the help of the perturbed paths and areas covered by VLF waves, the location of the earthquake was estimated.

Earthquake Temporal Prediction
Wavelet transform is used to detect any anomaly in radio data (Torrence and Compo, 1998). Earthquake temporal prediction (Wavelet analysis) has been used in many case studies related to radio wave perturbations and their relationship with earthquakes (Biagi et al., 2019;Righetti et al., 2012). Changes in VLF/LF amplitude data are very different from one path to another, so we need to analyse and compare data from different propagation paths in some ways. VLF/LF data should behave homogeneously when different propagation paths are analysed. To facilitate the comparison of different wavelet power spectra, it is desirable to perform a common normalization for the wavelet spectrum (Torrence and Compo, 1998). In this research, we have chosen the "Morlet function" as a wavelet. In this case, the wavelet transform of a time signal is a complex series that can be expressed by the square of its amplitude, i.e. we consider it as a Wavelet Power Spectrum. The implementation details are given in Torrence and Compo (1998).

Earthquake Spatial Prediction (Intersection-Union Method)
Before the occurrence of large earthquakes, VLF waves are perturbed in the transmitter-receiver paths. This perturbation occurs in paths where the earthquake is located in its fifth Fresnel zone and VLF waves are affected by ionospheric anomalies caused by the earthquake. Therefore, in an earthquake, not all of the paths are necessarily affected. Therefore, according to the area covered by all the paths and the separation of paths with and without perturbations, the probable area of the earthquake can be predicted. At this step, the overlay of existing layers using the proposed Intersection-Union method, the approximate location of the 2016 Kumamoto earthquake in Japan was determined. The implemented method in this research is a combination of intersection and union of areas based on the VLF/LF network geometry. To predict an earthquake by VLF waves, the earthquake must be within the fifth Fresnel zone of the receiver-transmitter paths. According to the undertaken studies in the Japanese VLF/LF network for more than 30 different earthquakes, if the distance of the epicenter to the receiver or transmitter is less than 300 km, the earthquake can be predicted (Hayakawa, 2015). Also, according to a research by Righetti et al. (2012) for the LF radio signals collected by the European VLF/LF network from 2009 to 2011, if an earthquake is within a radius of 300 km of the transmitter or receiver, it is possible to be predicted. Therefore, to estimate the location of the earthquake, the fifth Fresnel zone of all the receiver-transmitter propagation paths and the area around the receivers and transmitters within a radius of 300 km should be checked. In general, each receiver can receive signals from several transmitters, and each transmitter also sends signals to several receivers. Locations with a radius of 300 km around the receivers and transmitters are earthquake-prone areas only if there is a perturbation in all their propagation paths. In fact, in the fifth Fresnel zone of all the propagation paths related to that receiver or transmitter, perturbation should be observed; therefore, in a receiver or transmitter, even if one of the related propagation paths is not perturbed, a radius of 300 km is not considered and there would be no possibility of earthquake occurrence in that area. A VLF/LF network is a combination of the shapes presented in Figure (3). GIS provide the spatio-temporal capabilities that can be used in many disciplines (Bayat et al., 2020;Eslamimezhad and Delavar, 2019). In this research a method for approximating the probable location of an earthquake occurrence based on the general configuration of a VLF/LF network is proposed which is called Intersection-Union method. The steps of the proposed method are as follows: 1. Intersection of areas with a radius of 300 km around the receivers: For all receivers that receive signals from one transmitter and their received signal is perturbed, circles with a radius of 300 km around the receivers can only be assumed as the probable area of the earthquake occurrence that intersects to each other, and their overlapped area is considered as the area of the epicenter of the earthquake. 2. Intersection of areas with a radius of 300 km around the transmitters: For all transmitters whose signals are received by a receiver and the received signal is perturbed, circles with a radius of 300 km around the transmitters can be assumed as the probable area of the earthquake occurrence that intersects to each other, and their overlapped area is considered as the area of the epicenter of the earthquake. 3. Intersection of the fifth Fresnel zone of the propagation paths: The overlapped area among the different receivertransmitter paths is also the probable area of the earthquake occurrence and the intersection of these areas must be calculated. 4. The union of areas obtained from the steps 1 to 3 is the probable location of the epicenter. 5. The Intersection of the area obtained from the steps 1 to 3 is the probable location of the earthquake epicenter in the weighted form. Fresnel zones and 300 km areas around the receivers and transmitter in estimating the location of the earthquake epicenter have equal importance. The common parts of the circular 300 km areas around the receivers and transmitters has a weight of one, the common area of Fresnel zones has weight of one and the common area of Fresnel zones with circular 300 km areas has a weight of two.

Wavelet Analysis
According to Figure ( The results of the wavelet analysis are illustrated in Figures (5) and (6). The horizontal axis is the date (days) in 2016 at UT, the vertical axis is the Fourier period (minutes) and the colours of the image indicate the signal strength at different times. The time of the earthquake occurrence is indicated by an arrow.  Figure 6 shows the results of the wavelet analysis for shorter propagation paths over a period of more than one month (March 3 to April 18). In Figure 6 for the precursors with the shorter propagation paths, the anomaly is relatively obvious for all of these paths (JJI-ANA, JJI-IMZ, and JJI-TYH).
In Figure ( (6a) and (6b). The JJI-IMZ path showed a clear anomaly in early April 2016, as well as a clear anomaly from April 3 to 7, 2016 for the JJI-TYH. As observed, all of the propagation paths are perturbed and shown anomalies related to the earthquakes, but longer paths have more obvious anomalies. According to the results obtained to predict the time of an earthquake by the wavelet analysis method, the anomalies were observed about two weeks before the earthquake and one week before that had the highest amount of the anomaly.

The Intersection-Union Method
In this study, we have one transmitter and seven receivers. The transmitter sends signals to all of the seven receivers, so it has seven propagation paths. Each receiver also receives a signal from only one transmitter, so the receivers have only one propagation path. As previously mentioned, the perturbation has been observed in all the transmitter-receiver paths. Therefore, since there is no path for the transmitter or receivers without perturbation, all areas within a radius of 300 km around the transmitter and receivers and the fifth Fresnel zone are subject to earthquakes occurrence. Based on the network geometry, the steps for implementing the proposed Intersection-Union method are as follows: 1. Investigate and calculate the radius of 300 km (Righetti et al., 2012) around the receivers. 2. Investigate and calculate the intersection of the fifth Fresnel zone of the receiver-transmitter paths. 3. Determine the union of areas obtained at the steps 1 and 2 with an area of 300 km around the transmitter. 4. Weight the common areas in step 3. The circular 300 km area around the transmitter has a weight of one, the common area of Fresnel zones has weight of one and the common area of Fresnel zones with circular 300 km areas around the transmitter has a weight of two. In GIS, spatial data integration is usually undetaken with overlay analysis. Overlay analysis takes one or more datasets associated with a particular location as input and generates a new dataset as output (Yang et al., 2018). In GIS, Map Algebra functions such as polygon set operators (e.g., intersection and union) are applied to vector systems. The results of the implementation of the Intersection-Union method are shown in Figure (7). Figure (7a) shows areas with a radius of 300 km around the transmitter and receivers; Figure  (7b) illustrates the intersection areas with a radius of 300 km for the receivers and transmitters where Figure (7c), represents the area created by the intersection of the fifth Fresnel zone in Figure (2). Figure (7d) presents the union of intersecting the fifth Fresnel zone and a radius of 300 km around the transmitter, as the epicenter of the earthquake. As seen from Figure (7a), not all the circles related to the receivers have intersection in one common area, and the intersection of these seven areas is empty. Because all the circles of the receivers must have one area in common, not just two or more circles. For example, the circle of the STU and ANA receivers have nothing in common. Therefore, none of these areas can be a possible location for the epicenter of the earthquake. The result according to this method for the 2016 Kumamoto earthquake in Japan is an almost circular area with a radius of about 300 km as seen in Figure (7d). Due to the equal importance of Fresnel zones and 300 km areas around the receivers and transmitter in estimating the location of the earthquake epicenter, the same weight can be given to each one. A circular area with a radius of 300 km has weighed one, the common area of the fifth Fresnel zones has weighed one, and their common area has weighed two. Therefore, the yellow area in Figure (7d) (the common area of the fifth Fresnel zones) that is located inside the circle with radius of 300 km has a weight of two and it is the most probable area for the earthquake epicenter location. As can be seen in Figure (7), the epicenter is located very close to this area. Figure 7 (a) The areas with a radius of 300 km around the transmitter and receivers (b) The intersection of the areas with a radius of 300 km for the receivers and transmitter (c) The intersection of the fifth Fresnel zones (d) The union of intersection of the fifth Fresnel zones and a radius of 300 km around the transmitter that is the location of the epicenter of the earthquake. Hayakawa and Asano (2016) performed a preliminary analysis for the 2016 Kumamoto earthquake, and the VLF data were only studied for the JJI-IMZ path to investigate the presence of anomalies for the VLF emission before the earthquake. Asano and Hayakawa (2018) have also investigated the spatiotemporal evolution of lower ionospheric perturbations for the 2016 Kumamoto earthquakes by the VLF data observed in all the receiver stations. The results of the present study are consistent with those of Hayakawa and Asano (2016) and Asano and Hayakawa (2018).

Analysis of Temporal Results
Figure (8) illustrates the earthquakes that have been occurred between April 14 and 21, 2016 whose vertical axis present the intensity of the earthquakes. In addition to the main Kumamoto earthquake and its two pre-earthquakes, many earthquakes have been occurred. Although most of the earthquakes have a magnitude less than M= 6 and have no significant effects on VLF waves, several earthquakes between M=5.5 and M=6 may affect VLF waves in some propagation paths. Therefore, observing the difference in the time of observing anomalies in different paths in addition to the condition of the paths can also be caused by these earthquakes. Therefore, it would be difficult to claim with some certainty which day in each path corresponds to the observed precursor related to the earthquake. Therefore, according to the observed anomalies and previous statistical studies on the precursory with VLF waves, a period of about two weeks was obtained as a temporal interval for estimation of the time of earthquake occurrence.

Figure 8
The earthquakes that occurred on April 14-21, 2016 in the Kumamoto region of Japan (https://www.usgs.gov). (9) shows a map of the Kumamoto earthquake on April 15, 2016. As can be seen, the earthquake was felt with relatively high intensity in a large area with a radius of about 100 km. It is also felt with less intensity in an area with a radius of more than 300 km. Therefore, the estimation of an area with a radius of 300 km for the location of the epicenter of this earthquake has been done correctly. Also, the area with high intensity shaking has been overlaid by about 50% with the common area of the fifth Fresnel zones and the 300 km area around the transmitter. Therefore, the estimation of this area as the location of the epicenter with more weight has been correctly predicted. When there is only one receiver-transmitter path, the fifth Fresnel zone and areas with a radius of 300 km around the receiver and transmitter are the probable earthquake occurrence areas. For example, in this research, assuming the existence of only one path, the following results are obtained: The Shortest path (JJI-ANA): Assuming that only this shortest path exists, areas with a radius of 300 km around the receiver and transmitter and the fifth Fresnel zone of the receivertransmitter path, an area with about 1000 km length is the area prone to earthquakes. The Longest Path (JJI-STU): Assuming that only this longest path exists, areas with a radius of 300 km around the receiver and transmitter and the fifth Fresnel zone of the receivertransmitter path is the probable area for epicenter location, so an area of about 2000 km length is the area prone to epicenter location. Assuming these two scenarios, we examine the accuracy of the proposed model. The accuracy of the proposed model is given in Table 1. In this research, accuracy has used to express the radius of the epicenter area prediction and adaptation to its true radius and location. Therefore, accuracy is the radius of the area predicted for epicenter location. Accuracy expressed as a percentage is match between model output and reality. In other words, the ratio of the area predicted correctly by the method to the real shaking area.

Figure
Accuracy of the proposed model ~300 km Shortest path accuracy ~1000 km Longest path accuracy ~2000 km Accuracy of Shaking area prediction ~90% Accuracy of high intensity Shaking area prediction ~50% Table 1 The accuracy of the proposed model.
In this study, considering the data of a single transmitter and the network arrangement that all the receivers are located on one side of the transmitter, the accuracy of 300 km obtained is a good estimate of the earthquake epicenter location. There was also an earthquake that happened near the transmitter, affecting all of the paths. Therefore, it was impossible to further reduce the radius of the estimation area. The areas covered by the transmitter-receiver paths are very large, so the accuracy obtained is appropriate given the size of the area covered and the length of the paths. Short distances are more important in predicting the location of an earthquake. When the VLF method is used to predict the location of the earthquake, the entire perturbed area covered by a receiver-transmitter path is considered as a possible earthquake epicenter. Therefore, the area covered by the shorter paths is smaller for that a better spatial accuracy is achieved. As mentioned in the discussion of earthquake time prediction using the wavelet analysis, the longer paths showed clearer anomalies.
Therefore, the shorter paths are suitable for estimating the epicenter of an earthquake and longer paths are suitable for estimating the time of earthquake occurrence. In this study, by integrating the shortest and longest distance paths, both the location and time characteristics of the earthquake were estimated. To the best of the authors' knowledge, no research is reported so far on earthquake epicenter location analysis using VLF/LF data and method. Due to the significant reduction of the radius of the predicted area by the proposed method, the spatial uncertainty has decreased compared to other scenarios undertaken in this research. The impact of land use/land cover characteristics, topography of the area under study, the spatial topology of the transmitters and receivers, fusion of VLF/LF data with different earthquake precursors, implementing different spatial data fusion techniques for the earthquake epicentre location and time occurrence prediction, will be considered as further steps in the development of the proposed method.