Automatic hypocenter determination with the IPFx method for the 2018 Hualien earthquake sequence

The extended integrated particle filter (IPFx) method is an automatic source determination algorithm designed for the Japanese earthquake early warning (EEW) system. The method improved earthquake source determination during active seismicity by incorporating the smart phase association scheme. We applied this method to the 2018 Hualien earthquake sequence and evaluated its performance by comparing it to the manual catalog. We used 1-month continuous waveforms from February 2018 at 170 stations. Owing to the higher noise level, we improved the phase association algorithm to avoid noise contamination. Out of 127 earthquakes with a seismic intensity ≥ 4, 105 were successfully detected in one month, of which 103 had good accuracy with a location error of < 30 km. The detectability of earthquakes decreased immediately following large events. The IPFx method showed good performance in detecting earthquakes with seismic intensity ≥ 4 during the 2018 Hualien earthquake sequence. The method was also applied to the 1-day continuous data on April 18, 2021, and detected 14 earthquakes with a magnitude ~ 2 that were not on the manual catalog. Currently, the Central Weather Bureau in Taiwan uses the effective epicenter method to locate earthquakes for the EEW system. However, source determination for offshore events is difficult as most of the stations are on land. We expect the IPFx method to provide better location estimates for offshore earthquakes and during the period of active seismicity. It also provides an earlier warning as it sends the first message when three stations are triggered. This new method can potentially improve the speed and accuracy of the Taiwanese EEW system.


Introduction
The Central Weather Bureau (CWB) in Taiwan uses a regional earthquake early warning (EEW) system to warn the general public since 2016 . The system locates an earthquake based on the trigger information and estimates the magnitude from the displacement amplitude recorded at the near-source stations (Chen et al. 2019a, b). The seismic intensities (SI) at the target areas are estimated based on the ground motion prediction equation. To reduce the processing time, CWB uses the effective epicenter method (Chen et al. 2019a, b) to locate epicenters at a fixed focal depth of 10 km. The effective epicenter method estimates the location at the centroid of all triggered stations, and is effective owing to the simplicity of the calculations involved. For inland earthquakes, this method can provide an accurate estimation owing to the densely located seismic stations and good azimuthal coverage. However, as all the stations are inland, there is a large location error for offshore earthquakes (> 20 km, according to Hsu et al. (2018)) and thus large errors in the SI estimation.
Recently, two earthquake swarms occurred in the Hualien area on the east coast of Taiwan island. The 2018 Hualien earthquake sequence started on February 3. Beginning with an M 4.8 earthquake, 18 earthquakes of M 4.5 and larger occurred before the mainshock (the

Open Access
Terrestrial, Atmospheric and Oceanic Sciences *Correspondence: masumi@eqh.dpri.kyoto-u.ac.jp largest event in the sequence) on February 6 15:50 UTC (M 6.3) (Chen et al. 2019a, b). The active seismicity continued for a few days following the mainshock. On April 18, 2021, two earthquakes with M 5.8 and M 6.2 occurred at Shoufeng, 20 km southwest of Hualien. Active seismicity continued for approximately three hours.
Locating earthquakes during such active seismicity is challenging owing to phase association. Immediately following a large earthquake, multiple earthquakes occasionally occur in spatio-temporal proximity. In such cases, phases from two earthquakes co-occur, making it difficult to separate them. The new source determination method, the extended integrated particle filter (IPFx) method, incorporates the smart phase association scheme and works well for sequential earthquakes (Yamada et al. 2021). A similar method is currently used in the Japanese EEW system (Liu and Yamada, 2014;Tamaribuchi et al. 2014).
In this study, we applied the IPFx method to the continuous seismic waveforms during the 2018 Hualien earthquake sequence. We applied the method to the 1-month continuous waveforms in February 2018 and refined the parameters and phase association scheme. We tested the constructed IPFx method on 1-day continuous waveforms from April 18, 2021, and compared it with the CWB manual catalog for performance evaluation.

Data
We used; (1) 1-month continuous waveforms in February 2018 and (2) 1-day continuous waveforms on April 18 2021, recorded by the CWB (Central Weather Bureau Seismographic Network; CWBSN) and Academia Sinica, Institute of Earth Sciences (Broadband Array in Taiwan for Seismology; BATS), in total 170 stations with 3 component sensors as shown in Fig. 1. We used dataset (1) to tune the parameters in the algorithm and evaluate the performance and dataset (2) to test the fine-tuned algorithm. The CWBSN includes short-period, broad-band, and strong-motion sensors. If multiple sensors were installed at a station, we selected one sensor (the priority was a strong-motion, broad-band, or short-period sensor). To exclude events that occurred outside the target area, we removed the stations outside Taiwan (e.g., ocean bottom seismometer stations).
The CWB transmits the warning message to public via cell phone and TV. The public warning system (PWS) was developed and constructed by the government and communication company, and is based on the Cell Broadcast Service on the 4G network. The message is sent if the maximum estimated SI ≧4 and the estimated magnitude ≧5. The SI scale in Taiwan is a function of PGA, from zero to seven (Chen et al. 2019a, b).
The performance evaluation was carried out by comparing the CWB manual catalog and the IPFx catalog at convergence time. We selected the 127 earthquakes that had an observed SI of 4 or larger on the CWB scale in February 2018 (see Fig. 1a). We used the manual catalog instead of the CWB EEW results, as the EEW catalog is incomplete for the comparison. The result of the CWB EEW is only used for the time-history evaluation of the 2018 mainshock (Sect. 4.4). We focus on the earthquakes that occurred in the area around Taiwan (120 < longitude < 123 and 22 < latitude < 25). Most are shallow earthquakes with a focal depth of < 30 km.

A summary of the IPFx method
The IPFx method is a novel source determination algorithm developed for the EEW system (Yamada et al. 2021). The algorithm comprises two steps: single-station processing and network processing. The single-station processing step extracts station trigger and amplitude information from continuous waveforms (Yamada and Mori 2022), following which the network processing steps perform the source estimation from the extracted trigger and amplitude information based on Bayesian inference (Yamada et al. 2021). Both steps are performed in the central server so that we can adjust any types of data formats and change the parameters easily.
A special feature of the IPFx method is the smart phase association algorithm. During a period of active seismicity, such as that immediately following a large earthquake, phases from multiple earthquakes co-occur. This phase information has to be separated to accurately locate an earthquake. The phase association is performed based on the expected amplitude and arrival time.
When a single trigger is detected, we wait for another trigger near the first trigger station to start source determination. This status is referred to as a "pending earthquake". A predefined station group around the first trigger station is called a "trigger group". If the "virtual P-wave" passes the most distant station in the trigger group without additional triggers, the pending earthquake is deleted.
When multiple triggers (generally 3 triggers in the trigger group) are detected and the source determination process is started, any newly arrived trigger is evaluated to check whether it corresponds to the ongoing earthquake. This filter is very important: if noise (or a phase from another earthquake) joins the source determination process, the accuracy of the estimated source greatly decreases. We combined 4 parameters to create this filter: (1) the order of the trigger arrival, (2) the distance from the first trigger station, (3) the station (adjacent Voronoi cell) adjacent to the first trigger station, and (4) We refined this filter for use in this study (see the next section for details) as the noise level was higher than in the Japanese dataset in Yamada et al. (2021). The trigger information (P-wave arrival time and waveform amplitude) classified as this ongoing earthquake was used for source estimation (Yamada et al., 2021).

Improvement for the Taiwan network
We modified the original IPFx method designed for the Japanese network. As the noise level of the CWBSN dataset is higher than that of the Japanese network, the IPFx method detected many false events if we use the same parameters. Therefore, we found a stricter phase association algorithm through trial and errors. We reduced the number of stations in the trigger and estimation groups. The trigger group originally included all stations in the adjacent Voronoi cells as well as those within 30 km of the first trigger station. We reduced this distance from 30 to 20 km and the number of stations in the estimation group from 50 to 30. This resolved the noise contamination at the distant stations as these no longer contributed to the source estimation.
The original phase association algorithm is based solely on the error of the arrival time from the theoretical value. The theoretical arrival times are computed using the source location for the ongoing earthquake. If the error is less than 4σ (5 s, according to Yamada et al. (2021)), the trigger is included in the source estimation. We applied extra filters for the Taiwanese network as shown in Fig. 2. The fourth to sixth triggers must be either in the adjacent Voronoi cell or the closest 20 stations from the first trigger station with an arrival time error of less than 8σ (10 s). If more than six triggers are used for the source estimation, the location uncertainty is relatively small. Therefore, we applied a stricter condition for the seventh or later trigger arrival. The error of the station in the adjacent Voronoi cell and all other stations must be less than 8σ and 4σ, respectively, to be included in the source estimation.

Results
We detected 2618 earthquakes from the 1-month continuous waveforms in February 2018. Figure 3 shows the magnitude histograms of the IPFx and CWB manual catalogs. Earthquakes in the manual catalog with magnitude ≥ 5 were all detected. Earthquakes with a magnitude between 2 and 5 were missed if they occurred immediately after larger earthquakes. The detectability of  the earthquakes with a magnitude ≤ 2 was considerably lower. Although the sensitivity to small earthquakes may not be important for the purpose of issuing public warnings, it is essential for preventing false alarms.

Performance of the 2018 mainshock
We compared the performance of the IPFx method with that of the CWB EEW method. Figure 4 shows the result of the Hualien mainshock that occurred at 15:50:42 on February 6, 2018. The maximum observed SI was 7, and the CWB provided a public warning. The IPFx method does not consider the data transmission latency as it is an offline test. The earthquake occurred near the east coast of Taiwan. Although the station distribution is one-sided, the azimuthal coverage is reasonably good. The event was detected 4 s after the origin time, i.e., three stations were triggered within 4 s, and the maximum estimated SI exceeded the public warning threshold in 2 s. The public warning would have been issued 2 s after event detection. The CWB EEW system detected this event 13 s after the origin time and the public warning was provided 5 s after the detection. The reason for the delayed warning was that, at that time, the CWB EEW system required six triggers to alert the public. Both methods performed well at estimating location and magnitude. Figure 5 shows the result of the small aftershock (M 4.3) that occurred at 21:02:50 on February 6, 2018. The IPFx method determined the earthquake's location to be approximately 160 km away. The main reason for this inaccurate location was contamination by a minor earthquake not on the CWB manual catalog. Figure 6 shows the waveforms as a function of the distance of the M 4.3 earthquake. The closest station (HWA) did not operate at the time and the data was zero. The first trigger station (NACB) was likely noise, the No.2 trigger was from a minor earthquake, and the No.3 trigger was from the target earthquake. Therefore, the first estimation was based on a combination of the two earthquakes and noise. Triggers 2, 4, 9, 10, 12, and 13 were from the minor earthquake, and the others were from the target earthquake. In general, although the IPFx method was able to classify two events, these two events occurred too close spatio-temporally (within 2 s and 50 km). If HWA, the closest station, was operational, we would have been able to use the correct three triggers (HWA, ETM, and TWD) of the target earthquake instead of this wrong combination to start the source estimation process. Note that the magnitude of this earthquake was under the threshold (M = 5), so there was no public warning by the CWB EEW system. The number at the trigger line indicates the order in which the triggers contributed to the source estimation. Triggers Nos. 1, 2, and 3 were used for the event detection. The IPFx method scans triggers not used for the event detection for 3 s before the event detection, so Triggers Nos. 11, 12, and 13 were joined later. Triggers Nos. 5 and 7 were detected after the event detection and moved forward by trigger refinement

Performance of the 2018 Hualien earthquake sequence
We used the 127 earthquakes that had an observed SI ≥ 4 in February 2018 to evaluate the performance. The IPFx method detected 105 earthquakes, of which 103 had a location error of < 30 km. Figure 7 shows the time frame of the earthquakes detected by the IPFx method (bottom) and CWB manual catalog (top). Most of the earthquakes missed by the IPFx method occurred immediately following the large earthquakes on February 4 and 6. The details are discussed in the next section. Figure 8 shows the histogram with the location, magnitude, and origin time errors of the 105 detected earthquakes. Two earthquakes had a location error of ≥ 30 km, three earthquakes had a magnitude error of ≥ 1, and seven earthquakes had an origin time error of ≥ 4 s. The errors follow a Gaussian distribution with a mean of zero. There were six earthquakes where the IPFx method estimated the origin time approximately 10 s earlier than that of the CWB manual catalog. These cases were consecutive earthquakes that occurred at almost the same location. The first event was detected by the IPFx method; however, the second event was missed because the P-arrival was contaminated by the later phase of the first event. For EEW purposes, these differences in origin time are not critical as long as the magnitude and location estimations (and SI estimation) are reasonably accurate.

Performance of the 2021 Hualien earthquake sequence
We tested the developed method with the 1-day continuous data from April 18, 2021, which was the beginning of the 2021 Hualien earthquake sequence. There were two earthquakes exceeding the public warning threshold, and they were properly detected (open circles in the upper half of Fig. 9). To check the detectability of smaller earthquakes, Fig. 9 shows the time history of the earthquakes detected by the IPFx method and the CWB manual catalog. The manual catalog shows 294 earthquakes and the IPFx method detected 52 earthquakes. Most of the missed earthquakes occurred between 14:00 and 17:00, immediately following the two large events (M 5.8 at 14:11 and M 6.3 at 14:14). Fifteen earthquakes were detected by the IPFx method that were not on the CWB manual catalog. We checked the waveforms and determined that 14 were real earthquakes (see Additional file 1). Detectability of the IPFx method before 14:00 is in agreement with the CWB manual catalog.

Detection limit
As we have seen in Figs. 7 and 9, the detectability of earthquakes decreased immediately following a large event. We checked the relationship of the magnitude of the previous event and the time interval between the target and previous earthquakes. We expect the blind time (the period when no aftershock is detected) to be longer for an earthquake with a larger magnitude. Figure 10 shows the relationship of the time interval and magnitude of the previous earthquake for the 127 earthquakes with an SI ≥ 4. The blind time for the M 6.3 earthquake is approximately 7 min, whereas that of the M 4 earthquake is several seconds. Notably, the earthquakes can be detected even if the interval is less than 10 s, as long as the magnitude of the previous earthquake is small (M < 3). The earthquakes with an interval < 5 s are difficult to separate as any signals within 5 s after a trigger will be ignored by the phase detection algorithm (Yamada and Mori 2022). Note that this blind time depends on the size of the following event. If the second event is much larger than the first event, the blind time will be shorter. Figure 10 shows the upper boundary of the blind time.  Page 9 of 10 Yamada and Chen Terrestrial, Atmospheric and Oceanic Sciences (2022) 33:18

Magnitude difference
The magnitude of the IPFx method is estimated by the attenuation relationship of the JMA EEW (Liu and Yamada 2014). The magnitude is computed from the maximum amplitude of the vector sum of the displacement, hypocenter distance, and depth of the event. Therefore, the definition differs from that of the CWB catalog, which uses short-period seismograms (Shin et al. 2013). We compared the IPFx and CWB magnitudes for the 105 earthquakes in Fig. 11. The CWB magnitude is slightly larger for small earthquakes with M < 4.5; however, the two magnitudes agree reasonably well for large earthquakes. For simplicity, we used the IPFx magnitude for comparison with the manual catalog in this study.

Summary of the performance of the IPFx method
The IPFx method performed well in detecting earthquakes with SI ≥ 4 during the 2018 Hualien earthquake sequence. It detected 105 out of 127 earthquakes over the month, of which 103 had good accuracy with a location error < 30 km. The 24 missed events occurred within 7 min of large earthquakes. The detectability of earthquakes significantly decreased in the blind time immediately after a large earthquake. This is a limitation of the source determination approach. The retrospective test of the one-day dataset on April 18, 2021, demonstrated that the detectability of the IPFx method before the mainshock is in agreement with the CWB manual catalog. The IPFx method detected 14 earthquakes with a magnitude of about 2 that were not on the manual catalog.
The IPFx method detected the Hualien mainshock much faster than the current CWB EEW. The IPFx method would have been able to provide a warning 6 s after the origin time, whereas the CWB EEW required 18 s. The effective epicenter method locates an earthquake at the centroid of all triggered stations which leads to a large location error for offshore earthquakes (> 20 km, according to Hsu et al. (2018)). The average location error for the IPFx method is less than 10 km (Fig. 8a). Therefore, the IPFx method can, more accurately, locate offshore earthquakes where there are few stations and during the period of active seismicity when triggers from multiple events arrive simultaneously.

Remarks for real-time operation
Although this study is aimed for academic research, it is also meaningful to describe the operational system and possibility of future implementation. Currently, the data transmission from field stations to the central station is quite efficient in Taiwan. The data were digitized at a packet with a 1 s length and resolution of 100 samples per second. The data latency is approximately 1-2 s. All the real-time data streams were received by the software developed by Geotech, Guralp, and Nanometrics. To better integrate the data from different sources, the Earthworm software developed by the United States Geological Survey (Johnson et al. 1995) has been used, and the EEW system has also been developed on the platform.
The speed of the IPFx method is sufficient for realtime operation. The program can handle a dataset with approximately 1000 stations × three channels in real time. The computation time required for reading 1-h data from the hard disk and applying the IPFx method was 36 min (CPU: Xeon 3.46 GHz, OS: Linux CentOS 6, memory: 48 GB). To implement the IPFx program, we need to pre-process the CWB real-time data streams and extract trigger information. Currently, the real-time data streams are all integrated into the Earthworm platform. The challenge for the real-time operation is that the Earthworm software is difficult to modify because the source code was created using the C language. However, recently, a python interface to the Earthworm Messaging Transport system (PyEarthworm) was created (Hernandez et al. 2018). The Python language and its toolbox for seismology (ObsPy) simplify communications with the Earthworm system (Beyreuther et al. 2010). Therefore, they will increase the feasibility of improving the CWB EEW system with the IPFx method.

Conclusion
We improved the IPFx method used for the Japanese EEW system, applied it to the 1-month continuous waveforms of the 2018 Hualien earthquake sequence, and evaluated its performance by comparing it to the manual Fig. 11 Magnitude of the CWB manual catalog and IPFx method for the 127 earthquakes with an observed SI ≥ 4