High‐Resolution Fault‐Rupture Imaging by Combining a Backprojection Method With Binarized MUSIC Spectral Image Calculation

Backprojection (BP) methods are widely used for estimating the fault rupture processes; however, they are inherently susceptible to noise. Hence, noise suppression is an important research target. In this paper, we develop a fault rupture imaging method by combining beamforming‐based BP and MUltiple Signal Classification (MUSIC), which realizes artificial noise suppression at a high spatial resolution. The stations are grouped into arrays according to the SH wave radiation coefficients, and MUSIC analysis is performed on each array. The MUSIC spectral images of these arrays are binarized and then multiplied by the BP images. Spatial filtering is also applied to the images based on the possible range of the rupture velocity and rise time. When tested using synthetic test data, the proposed method worked as expected. We then applied this method to the 2016 Kumamoto earthquake by interpolating the travel times from the observed travel times of relocated hypocenters using a 3‐D velocity structure model. In the area of large slip and slip rate approximately 30–50 km from the hypocenter on the Futagawa fault, the spatiotemporal evolution of the fault ruptures and waveform inversion results were generally in harmony. The distributions of the low‐ and high‐frequency seismic radiations are complementary, as is understood in the context of fault rupture physics. This method can aid in understanding and modeling the details of seismic radiation sources, enabling the accurate prediction of strong ground motion even in near‐fault areas.

The BP method and its derivatives are now widely used as one of the most popular source process analysis methods. Although there are a variety of BP methods, they all share the same basic principle, which is essentially that of migration in the reflection method. By stacking the seismic records obtained from all of the stations in a region, the amplitude of the stacked seismic records increases at the location and time of actual seismic wave emittance. This method has been used in many studies because it is simple in principle, requires negligible a priori information for analysis, and is computationally inexpensive. Additionally, compared with waveform inversion, the fault rupture process can be examined more directly and at a higher resolution because the use of regularization in waveform inversion, such as spatiotemporal smoothing, is not essential for BP. Therefore, BP is used in various manners in a wide range of disciplines, such as investigations of the spatiotemporal evolution of a fault rupture and the immediate determination of the location of an earthquake occurrence.
BP methods are inherently susceptible to noise, such as artificial noise, known as "swimming effects" (Meng, Ampuero, Luo, et al., 2012). Various studies have attempted to obtain a high-resolution seismic source image while suppressing noise. Some of these methods include correction of the travel time for the seismic source location  or for the seismic source area using Kriging (Palo et al., 2014), travel time correction with aftershock data (Meng et al., 2016), the use of the travel time relative to the hypocenter (Takenaka et al., 2009;Yamamoto & Takenaka, 2006), the use of the travel time from a reference station (Maercklin et al., 2012), and the use of the travel time calculated with a 3-D velocity structure model (Liu et al., 2017;Song & Ge, 2019).
The grouping, selection, and weighting of stations are also important for improving rupture images; for example, grouping of stations according to their azimuth or the use of multiple arrays by grouping adjacent stations (Kiser & Ishii, 2012), selection of stations by semblance (Rössler et al., 2010), desampling of stations (D' Amico et al., 2010), selection of stations to minimize artifacts using genetic algorithm (Kehoe et al., 2019), weighting of stations according to the epicentral distance (e.g., Pulido et al., 2008), and weighting of stations according to their azimuthal distribution to the hypocenter (Jakka et al., 2010).
Based on these studies, the suppression of noise, in addition to improvements to the resolution of the BP methods, are important research targets. Furthermore, for the BP method using local strong-motion waveforms, the recorded seismic waves propagate through the strongly heterogeneous shallow part of Earth; thus, various 3 of 36 innovations are necessary to accurately image the fault rupture process. In this paper, we develop an imaging method for fault rupture by combining a BP method based on beamforming and MUSIC with multiple arrays. We test the proposed method by using a synthetic dataset. We then apply the method to the 2016 Kumamoto earthquake to confirm the effectiveness of the proposed method for practical data. We exploit realistic travel time data obtained by interpolations of the observed travel times that were calculated from the relocated events surrounding the causal faults of the 2016 Kumamoto earthquake, with the 3-D velocity structure reported in Matsubara et al. (2017). Figure 1 shows an overview of the proposed methodology. The stations surrounding the source fault are first grouped into multiple arrays based on the radiation coefficient of the SH wave at each station. SH waves do not interact with P or SV waves and thus are suitable for extracting information about the seismic source. By selecting stations with large SH wave radiation coefficients and grouping stations according to the sign of the SH wave radiation coefficients, we obtain information regarding the seismic source based on SH waves, which are relatively less distorted. Next, the 2-D MUSIC spectra are calculated along the fault plane (MUSIC spectral image) for each array. The MUSIC spectra are binarized, multiplied, and filtered based on the possible range of rupture velocity and assumed maximum rise time. Then, the BP image is obtained based on all the arrays used in the calculation of the MUSIC spectral images. Finally, the binary and BP images are conjoined. Using this MUSIC spectral image conjunction, we obtain the location of the seismic wave radiation source on the fault plane by suppressing the artificial noise at a high resolution.

Summary of the Proposed Method
MUSIC (Schmidt, 1986) can estimate the incidence angles of incoming waves at a high resolution; however, the accuracy of the amplitude estimation is low. In contrast, beamforming has a low resolution when estimating the direction of incident wave arrival, but has high accuracy when estimating the incidence wave amplitude (e.g., Kikuma, 2003). By combining MUSIC and beamforming, the advantages of each method can compensate for their respective disadvantages.
MUSIC uses the orthogonality between the eigenvector corresponding to the noise eigenvalue of the waveform correlation matrix inside an array and the array response vector. This subspace method produces high-resolution results. When a wave radiation source exists on a fault plane, a MUSIC spectral image with an array has peak values along the straight line connecting the array with the seismic radiation area on the fault plane. When a MUSIC spectral image is calculated for each of the multiple arrays, the direction of each array, as observed from the seismic radiation area, is different; thus, linear structures with large amplitudes in the MUSIC spectral image of each array intersect in the seismic radiation area on the fault plane. Therefore, the seismic radiation area on the fault plane can be detected by extracting the area where straight lines, with large amplitudes, in the MUSIC spectral images of all arrays intersect; within this area, the amplitude of the MUSIC spectrum of each array has a large common value. By multiplying the result of the beamforming-based BP by a weight of 1 for the area extracted via MUSIC and 0 for other areas, the quantity of seismic wave radiation along the fault plane can be obtained at a high resolution. Figure 2 illustrates the fundamentals of the proposed method.
One of the most significant drawbacks of the BP method is artifacts, and the proposed method can suppress them. In particular, as the swimming effect sweeps from the source to the station (e.g., Meng, Ampuero, Luo, et al., 2012) which varies from array to array, we can mitigate this noise by multiplying the MUSIC spectral images of multiple arrays. In addition to the swimming effects, artifacts caused by the local ground structure can also be suppressed using a multi-array. Seismic signals related to fault ruptures are common in seismic waveforms at all stations, while seismic noise, such as noise caused by the local ground structure, is not expected to be present in all the arrays. This seismic noise can be suppressed by multiplying the rupture images from all the arrays. The procedures of the proposed method are detailed in the following sections.

Calculation of MUSIC Spectral Images
To obtain a MUSIC spectral image, we perform eigenvalue decomposition of the correlation matrix among the waveforms of the seismometers that are composed of the array. The obtained eigenvalues are then divided into those corresponding to eigenvectors that span the signal subspace and those corresponding to eigenvectors that span the noise subspace. In other words, the number of incident waves in the array is determined. Then, the MUSIC spectrum is calculated from the array response vector and eigenvectors that span the noise subspace. To divide the eigenvalues into those corresponding to the signal and noise subspaces, the following methods are generally used: (a) the method that counts the number of eigenvalues normalized by the largest eigenvalue larger than a certain threshold (Asano, 2011), or (b) the method using Akaike's Information Criterion (AIC: Akaike, 1974), or (c) the method using Minimum Description Length (MDL: Rissanen, 1978). While MUSIC assumes the whiteness of noise, seismic motion consists of a complex wavefield containing noise correlated with signals. Teleseismic waveforms mainly consist of seismic phases that have propagated through the relatively The MUSIC spectral image on the fault plane has peaks on the line connecting each array and seismic radiation area on the fault plane. MUSIC spectral images are calculated for each of the multi-arrays. The seismic radiation area is extracted from the intersection of the large MUSIC spectral amplitude line from each array. The weights are set to 1 in the extracted seismic radiation area and to 0 in other areas (binarization). The right panel shows the beamforming step. Beamforming is executed using all the stations for all of the arrays, and the obtained rupture image is multiplied by the binarized weights obtained in the MUSIC step. The black squares represent stations that comprise each array. The blue and red colors on the fault plane signify the seismic radiation amplitude evaluated by beamforming.
10.1029/2022JB024003 5 of 36 homogeneous deeper part of the Earth, where the underground structure has a relatively minimal influence. In contrast, local strong-motion records include seismic waves which have traveled through the strongly heterogeneous shallow subsurface such that the underground structure significantly influences the wavefield. This results in complex wavefields with noise correlated to the signal components. Under such a situation, MUSIC does not function properly and the method that divides the eigenvalues also fails.
Eigenvalue decomposition is used in both MUSIC and principal component analysis (PCA: Pearson, 1901). One of the objectives of PCA is dimension reduction. In dimension reduction, eigenvalues are selected from among the eigenvalues obtained by decomposing the variance-covariance matrices, such that the data can be represented by a minimally optimal number of eigenvalues and their related eigenvectors without losing information present in the original data. In the PCA used to reconstruct data or perform dimension reduction, eigenvalues with a value of >1 or eigenvalues whose cumulative contribution ratio reaches a certain value, such as 80%, are selected (e.g., Asano, 2011). In this paper, the latter approach is adopted. Such an approach has already been applied in seismology; for example, to detect low-frequency tremor sources (Ueno et al., 2010). Eigenvalues corresponding to the signal subspace in the MUSIC spectral image are selected so that the cumulative contribution rate is 70%. In other words, the eigenvalues are selected in order from the largest to smallest until we obtain 70% of the sum of all of the eigenvalues. The cumulative contribution ratio is set to 70% because we found that 60%-80% gave reasonable results based on the preliminary tests.
Here, we present a synthetic example (Figure 3). We consider a fault and four adjacent seismic arrays. The fault plane is 22 km long, 12 km wide, and dipped at 45° in the positive direction of the y-axis. The four arrays have a minimum distance of 50 km from the center of the fault. Each of the four arrays consists of nine stations spaced 20 km apart (in accordance with the spacing used by the strong-motion seismograph network in Japan, https:// www.kyoshin.bosai.go.jp/kyoshin/) in both the direction of the x-and y-axis. There are three seismic radiation sources on the fault, indicated by stars, which successively radiate pulses at 1.0 Hz for 1.0 s. The time interval between the pulses emitted by the three sources, sources 1, 2, and 3 is 3.0 s. The wave propagation speed is assumed to be 3.0 km/s. Figure 4 shows an example of a waveform observed at a station. For each of the four arrays shown in Figure 3, the MUSIC spectrum is calculated on a grid along the fault plane to obtain a MUSIC spectrum image. In this case, the reference point used in MUSIC is source 1. In many cases, the hypocenter of an earthquake is determined immediately after the occurrence of the earthquake. Using this a priori information, we set the hypocenter location as the reference point for the MUSIC analysis.
Figures 5-7 show the MUSIC spectral images on the fault plane obtained for each array at 0.2 s after the start of emission from sources 1, 2, and 3, respectively. The MUSIC spectrum image has a peak on the line connecting the array and seismic source. Since the azimuth of the seismic source observed from an array varies from array to array, the location of the seismic source can be extracted via the conjunction of the MUSIC spectrum images obtained from each array and by extracting the peak common to the MUSIC spectrum images from all of the arrays.

Binarization of MUSIC Spectral Images
Since the amplitude of the MUSIC spectrum image varies by orders of magnitude, these images calculated for each array are first binarized and then multiplied to retain signal components with small amplitudes. If the MUSIC spectral images of all arrays are multiplied and binarized, small signal components in one or more of the arrays disappear if the spectral amplitude is small in one of those arrays. Under such a situation, we cannot obtain signal components with small amplitudes by setting a certain threshold and selecting values larger than that threshold. The amplitudes of the MUSIC spectral  images at sources 1-3, as shown in Figures 5-7, are not large compared with the maximum value on the fault plane. Therefore, we adopt the cumulative contribution concept in PCA: we sort the grids on the fault plane according to the MUSIC spectral amplitude and select amplitude grids from the largest to smallest until the sum of their amplitudes exceeded 70% of the total MUSIC spectrum image amplitude. We assign a value of 1 to grids corresponding to the selected grids and 0 to the remaining grids. This binarization operation also allows the identification of weak signal components and the detection of seismic wave radiation sources.

of 36
The left part of Figure 6 shows the binarized MUSIC spectral images on the fault plane from each array at 0.2 s after the initiation of emissions from sources 1, 2, and 3. The binarization operation detects the region containing the radiation source. By taking the product of the binarized MUSIC spectral images of each array, which is a type of logical conjunction, we can identify the radiation sources commonly captured in the MUSIC spectral images of all of the arrays while suppressing the artificial noise. The right part of Figure 6 displays the product of the binarized MUSIC spectral images on the fault plane obtained by each array at 0.2 s after the start of emission from sources 1, 2, and 3. Figure 6 illustrates that the region including the radiation source can be extracted by taking the product of the binarized images from all of the arrays.

Filtering Based on Rupture Velocity and Rise Time
To suppress the artificial noises further, spatial filtering based on the possible range of rupture velocity and assumed maximum rise time is performed as follows: where ( ) is the coefficient used to multiply the i-th grid at time t, and are the minimum and maximum values of the assumed rupture velocity, respectively, is the maximum value of the possible rise time, and d is the distance from the hypocenter to the i-th grid. By using the largest possible values for both and , as well as the smallest possible values, artifacts can be suppressed by avoiding the distortion of the rupture image by such a filtering process. The value of can be set with reference to the slip velocity time function inferred from waveform inversion, if available. If such information is not available in advance, the approximate value can be set based on the scaling relation between rupture time (fault length divided by rupture velocity) and rise time given by Kanamori and Anderson (1975). The rupture velocity often exceeds the S-wave velocity (e.g., Bao et al., 2019;Bouchon & Vallée, 2003;Dunham & Archuleta, 2004;Hicks et al., 2020;Kehoe & Kiser, 2020;Sangha et al., 2017;Walker & Shearer, 2009;Wang et al., 2012;Yue et al., 2013); however, for earthquakes not characterized by a supershear rupture, the possible range of rupture velocity values is usually capped at the Rayleigh wave velocity (Broberg, 1996). We exploit this prior information to suppress artificial noise and improve the source imaging accuracy for earthquakes without a supershear rupture. Figure 7 shows the results of spatial filtering based on this rupture velocity. In the example shown here, and are 1.0 km/s and 3.0 km/s, respectively, while assuming = 10.0 s. We can see that the accuracy of the seismic radiation estimation can be improved by imposing a constraint on the seismic source location in space and time based on the possible range of rupture velocity and rise time.

Conjunction of Binarized and BP Images
Although MUSIC has an excellent resolution for estimating the incident direction of arriving waves, it has poor accuracy for evaluation of the power of incident waves under a noisy environment; therefore, the latter must be estimated via another method that can accurately evaluate the power, such as beamforming (e.g., Kikuma, 2003). Many BP methods are based on beamforming. By combining MUSIC, which has superior resolution for estimating the direction of wave arrival, and BP based on beamforming, which has superior accuracy for evaluation of the power. The disadvantages of both methods can be compensated by the advantages of each method. We utilize all four arrays shown in Figure 3 to estimate radiation on the fault plane via BP based on beamforming.
The travel time used for BP is calculated with a wave propagation velocity of 3.0 km/s; the station correction is performed such that the first pulses from all of the stations are synchronized when stacked at source 1. Figure 6 shows that although all of the seismic sources (i.e., sources 1, 2, and 3) are detected, artificial noise is present surrounding the actual radiation sources, as indicated by the stars. In Figure 7, the left part shows the multiplied image (from Figure 6), which gives the final image after being subjected to the rupture velocity filter followed by multiplication with the BP image. This shows that using the information on the direction of incident wave arrival from MUSIC, the rupture propagation velocity, and rise time can suppress the artificial noise and improve the accuracy of source imaging by BP.

Tests With Noisy Data
In earthquakes, rupture velocity is not constant over the fault plane, and the seismic waveform records contain noise. In order to verify the effectiveness of the proposed method for such a realistic case, we use the same station deployment and fault setting as in the previous test shown in Figure 3, vary the rupture velocity on the fault plane, and add uniform random noise to the input waveform. The noise test is conducted for three cases with signal-to-noise ratios of −3.01 dB, −6.02 dB, and −10.0 dB. We show the results are shown on the left, center, and right of Figure 8, respectively. The uppermost panels of Figure 8 indicate the waveform data including noise employed in the tests for each case with signal-to-noise ratios of −3.01dB, −6.02dB, and −10.0dB, from the left to the right, respectively. In all cases, the rupture velocity has been set to 2.0 km/s from source 1 to source 2 and to 2.5 km/s from source 2 to source 3. We have set , , and , to 2.0 km/s, 3.0 km/s, and 3.0 s, respectively. Figure 8 illustrates that even for the case where rupture velocity changes on the fault plane, the radiation sources can be detected at T = 0.2 s, T = 4.3 s, and T = 7.5 s, which are the time source 1, source 2, and source 3 emitted the waves, respectively. We also see that even in the case with a signal-to-noise ratio of −6.02 dB, the present method identifies radiation sources at every radiation time.
The reason why the present method works well with noisy waveform data may be as follows: when the waveforms contain uniform random noise uncorrelated among the stations in the array, it only adds a constant value to the eigenvalues of the correlation matrix of the waveforms in the array and does not affect the eigenvectors. Thus MUSIC's performance is not degraded (Kikuma, 2003). For actual earthquakes, especially local events, seismic waves propagate in the shallow underground with large heterogeneity, resulting in a variety of later phases. These later phases can become noise that correlates with the signal, coming in from different directions than the true source. These noises arrive at the stations in the array shortly after the arrival of the true source signal, which will deteriorate the performance of MUSIC.

Tests With Radiation Pattern and Theoretical Waveforms
In Subsections 2.5 and 2.6, the underground structure was assumed to be homogeneous, and simple sine waves were employed without considering the radiation coefficient. However, the real underground structure is not homogeneous, the waveforms at the stations are not simple sine waves. Also, in a real earthquake, the radiation coefficient of the seismic wave at each station varies with the station azimuth. To test the present method in a more realistic setting, a 1-D velocity structure model is set up from the 3-D velocity structure model of Kumamoto by Koketsu et al. (2008) to calculate the synthetic waveforms, which include the azimuth-dependent 10.1029/2022JB024003 9 of 36 radiation coefficient at each station of three arrays as we will use for the analysis of the 2016 Kumamoto earthquake in the next section. Figure 9 shows the array and fault configurations; three arrays consisting of nine stations with 30 km intervals have been set up for testing. The strike, dip, and rake of the faults are set to 270°, 45°, and 180°, respectively. The synthetic seismograms are calculated using the reflectivity method based on the discrete wavenumber method of Bouchon (1981) and the reflection-transmission matrices of Kennett andKerry, 1979 (e.g., Takenaka &Sasatani, 2000;Takenaka & Watanabe, 2021;Takeo, 1985;Yao & Harkrider, 1983). In the calculation, a bell-shaped function which has the pulse width corresponding to each frequency was employed as the source time function. Figure 10 displays the waveforms at each frequency for the three arrays. As in Subsections 2.5 and 2.6, three seismic emission sources are placed on the fault plane: source 1, source 2, and source 3. The rupture velocity varies on the fault plane, propagating at 2.2 km/s from source 1 to source 2 and 2.6 km/s from source 2 to source 3. Here, tests are conducted at five frequencies: 0.25 Hz, 0. The left, center, and right column depicts the results for the waveform of the signal-to-noise ratio of −3.01 dB, −6.02 dB, and −10.0dB, respectively. The rupture velocity is 2.0 km/s from source 1 to source 2 and it changes to 2.5 km/s from source 2 to source 3. The amplitude has been normalized to the maximum value. The second upper panels display the image at time T = 0.2 s, where the wave was radiated by source 1. The third upper panels show the image of source 2 at time T = 4.3 s. The bottom panels show the image of source 3 at time T = 7.5 s, where the wave was emitted by source 3. The stars on the plane represent source 1, source 2, and source 3, from the left to right, respectively. and 2.6 km/s, respectively; the cases with and set to 1.0 km/s, 2.6 km/s, and 2.2 km/s and 3.0 km/s, respectively.
The results at each frequency are shown in Figure 11. The locations of the seismic emission sources are successfully imaged at each frequency although the high-frequency images have wide weak scattered ghosts for sources 2 and 3. The three seismic sources are captured at the correct timing, which indicates that the present method is also effective for imaging ruptures with varying propagation velocities on the fault plane. The size of the imaged radiation source depends on the frequency; at lower frequencies, such as 0.25 Hz, the size of the imaged radiation source is larger than at higher frequencies, such as 2.0 Hz, because the seismic wave wavelength becomes longer as frequency decreases if seismic wave velocity is the same. In Figures 7 and 8, the seismic emission source is clearly captured, whereas in Figure 11, the captured image is blurred and false images occur. This is because the radiation coefficient varies depending on the azimuth of the station from the source, and the underground structure generates various waves in addition to direct waves, which degrades the performance of beamforming and MUSIC. In particular, the generation of later phases that correlate with the direct waves and are mixed into the waveform will reduce the accuracy of MUSIC's excellent estimation of the direction of the incident waves. However, although false images with relatively weak radiation strength are scattered over the fault plane, the radiation strength is stronger (dark red) at the location of the seismic emission source, indicating the validity of the present method for practical application.
The proposed method images the rupture that propagates with an average velocity between e and in Equation 1. At each point on the fault plane, imaging begins at the time calculated by and ends at the time calculated by plus . Therefore, neither rupture propagating faster than nor rupture propagating slower than can be captured. The difference between the time defined by and the time defined by is longer at points farther away from the rupture initiation point, so false images are also more likely to occur. If is greater than the actual rupture velocity, the time when the rupture has not yet actually arrived will be imaged, resulting in a false image. Similarly, if is smaller than the actual rupture velocity, the time when the fracture has already been completed will be imaged, resulting in artifacts. In the case where and are set to 2.6 km/s, the difference from the actual rupture velocity is small, so the length of time covered by imaging is not too long or too short at each point on the fault plane, and the artificial noise is suppressed. On the other hand, in the cases where and are set to 1.0 km/s and 2.6 km/s, respectively, is smaller than the actual value, and the time length of imaging at each point on the fault plane is too long, resulting in more ghosting. For the cases where and are set to 2.2 km/s and 3.0 km/s, respectively, the values are larger than the actual ones, resulting in too long a time for imaging at each point on the fault plane, and thus more artifacts. These results indicate that although artifacts occur when and differ significantly from the actual average rupture velocity, the seismic emission sources are approximately successfully captured in the artifacts distributed on the fault plane as images with high radiation strength. Therefore, as long as and do not deviate too much from the actual average rupture velocity, the proposed method is robust to the values of and .
The closer both and are to the actual rupture velocity, the less false images are produced. The value of can be given approximately by the empirical relationship between rupture time and rise time of the fault by Kanamori and Anderson (1975), as described before. If is smaller than the actual rupture velocity, or if the value of vrmax is larger than the actual rupture velocity, the time that each point on the fault plane is subject to imaging is longer than the time it is radiating seismic waves, resulting in many false images. On the other hand, if is larger than the actual rupture propagation velocity or is smaller than the actual rupture propagation velocity, the time when each point on the fault plane is emitting seismic waves is not subject to imaging, and the seismic wave radiation source is not captured correctly. Thus, since the results from the proposed method are controlled by the settings of and , the accuracy of the obtained results can be improved by making the values of and as close to the actual rupture velocity as possible. For this purpose, it is recommended to use a priori information such as the triggering speed of the first-time window in the waveform inversion, if  available. If such a priori information is not available, the average S-wave velocity in the source region could be used with an empirical relationship between S-wave velocity and rupture velocity (e.g., Geller, 1976).

2016 Kumamoto Earthquake
The 2016 Kumamoto earthquake occurred at 1:25:5.4 (Japan Standard Time: +9 h to Coordinated Universal Time) on 16 April 2016, along the Futagawa-Hinagu fault system, which is part of the Median Tectonic Line, a large fault system that runs through the Japanese archipelago. According to National Research Institute for Earth Science and Disaster Resilience (NIED), the focal mechanism was a right-lateral strike-slip, and the Japan Meteorological Agency (JMA) magnitude (Mj) was 7.3. Two days before this earthquake, an earthquake of Mj 6.5 occurred while an earthquake of Mj 6.4 occurred the day before. Including the Mj 7.3 event, the Futagawa-Hinagu fault zone experienced eight earthquakes of ≥Mj 5.8 and seven events with JMA seismic intensity scale of 6 lower, which struck the surrounding area over five days from April 14 to April 18. The 2016 Kumamoto earthquake and its series of seismic activities have yielded an awareness of the danger of repeated strong ground motion events over a short period.
The 2016 Kumamoto earthquake occurred in the volcanic zone of the Aso volcano, which is an active volcano; the Futagawa fault was recognized as it does not traverse the area inside Mt. Aso. However, displacement along the Futagawa fault traversed through the caldera of Mt. Aso during the 2016 Kumamoto earthquake (Ozawa et al., 2016). This earthquake showed that fault plane setting for future earthquakes in volcanic regions must take into account the fact that faults can move and cause large seismic motions even inside volcanic bodies. The occurrence of this earthquake has also demonstrated that even active faults inside a volcano can cause large earthquakes.
Diverse studies have been conducted on the 2016 Kumamoto earthquake, such as a survey of the surface displacement (e.g., Toda et al., 2016), aftershock distribution (e.g., Kato, Fukuda, et al., 2016;Uchide et al., 2016;Yano & Matsubara, 2017), crustal deformation analysis, such as InSAR analysis and pixel-offset analysis of SAR images(e.g., Himematsu & Furuya, 2016;Ozawa et al., 2016), teleseismic waveform inversions (e.g., Yagi et al., 2016), and inversions of strong ground motion records (e.g., Asano & Iwata, 2021;Kobayashi et al., 2017;Kubo et al., 2016;Uchide et al., 2016), among others. Slip partitioning has been reported for the short faults parallel to the Futagawa fault, where the main slip occurred (Himematsu & Furuya, 2016;Toda et al., 2016). The 2016 Kumamoto earthquake was a complex earthquake with multiple source faults. Since the 2016 Kumamoto earthquake occurred in a well-developed strong-motion network, we obtained abundant strong-motion records in both the azimuth and distance distributions, which are appropriate for verifying the effectiveness of the proposed method. Therefore, we estimated the fault rupture process of the 2016 Kumamoto earthquake with the proposed method using strong-motion records.

Array Configuration
In the proposed method, multiple arrays must be placed around the source fault. The 2016 Kumamoto earthquake is an inland event, and its causal faults are all surrounded by stations. The focal mechanism is mainly of a right-lateral displacement type, where the area surrounding the causal faults could be divided into four regions according to the radiation coefficients of the SH waves. Therefore, it should be possible to construct four arrays by grouping together the stations that exist in each region. However, in the northeastern part of the epicenter of the 2016 Kumamoto earthquake (Oita area), a sufficient number of stations are not available for our analysis; hence, the selected stations were divided into only three groups (FKOH, MYZH, and KGSH), and three arrays were constructed. In the Oita area, 80 km from the 2016 Kumamoto earthquake, the strong ground motion radiating from the 2016 Kumamoto earthquake triggered an M6-class event (Nakamura & Aoi, 2017;Yoshida, 2016). The M6-class event occurred 33.5 s after the origin time of the 2016 Kumamoto earthquake  (Nakamura & Aoi, 2017) and that contaminated the waveform data in the Oita area. Thus, we excluded stations in the Oita area. Figure 12 shows the configuration of the three arrays.
In general, the greater the spatial extent of the array, the better the spatial resolution, and the greater the number of seismometers in the array, i.e., the higher the spatial density of stations, the better the signal-to-noise ratio. On the other hand, to accurately estimate the fault rupture process using BP based on beamforming and MUSIC, the seismic waveforms at the stations composing the array must retain a common signal that contains information on the fault rupture process at a high signal-to-noise ratio. In other words, the seismic waveforms of multiple stations that compose the array must be highly correlated. If the spatial extent of the array is too large, the subsurface structure and seismic wave propagation paths will differ significantly among the stations in the array. As a result, the waveforms of the stations in the array will differ greatly from each other, causing destructive interference in waveform stacking in the beamforming and correlation matrix calculation in MUSIC, making it impossible to accurately extract information on the fault rupture process. Therefore, it is necessary to set the spatial extent of the array so that the site and path effects are similar for each station in the array, resulting in a high correlation among waveforms recorded by the stations.
In addition, for seismic waves emitted from any point on the fault plane, the waveforms observed in the array must have the same polarity, which also imposes restrictions on array configurations. Keeping this in mind, stations close to the nodes of fault planes were excluded, which was done by setting a threshold of 0.1 for the radiation coefficient and removing stations with radiation coefficients smaller than the threshold. This is because the observed polarity near the nodes of radiation may differ from the calculated value owing to the complexity of the earthquake mechanism, estimation errors associated with the mechanism, and the influence of the local subsurface structure near the stations. If the actual polarity is different from the calculated value, destructive summation occurs among the waveforms that compose the array, and the signal component that contains information on the fault rupture is lost.
In addition to these conditions, the stations in the array must be unified, either near the fault where the near-field term dominates or far from the fault where the far-field term dominates. This is because the shapes of the waveforms of the same earthquake, in which the near-field term predominates, and in which the far-field term predominates differ significantly. In general, it is rare to have enough stations in the immediate vicinity of a fault where the near-field term dominates. In addition, since the far-field term in velocity records is proportional to the slip acceleration of the fault, the physical interpretation of BP results obtained by stacking seismic waveforms in which the far-field term dominates may be easier, so stations in the far-field may often be employed. The present method uses MUSIC, and in MUSIC, the maximum number of identifiable incident waves is the number of stations in the array -1 (e.g., Kikuma, 2003). Hence the minimum number of stations comprising the array should be the expected number of incident waves to the array +1 or more. MUSIC takes the harmonic mean of the angular spectra, the number of which equals that of stations in an array minus that of incident waves. As a result, signal components common to all angular spectra are extracted, false images are suppressed, and stable and accurate direction estimation is achieved (e.g., Kikuma, 2003). In general, the more stations than the number of incident waves +1, the more MUSIC benefits are exploited and the better the performance. Furthermore, in practice, not Figure 11. Results for waveforms calculated for 1-D velocity structure with considering radiation coefficients. Results at 0.25 Hz, 0.5 Hz, 1 Hz, 2 Hz, and 4 Hz are illustrated from left to right. They are plotted at the time when the amplitude of each frequency wave emitted is at its maximum, i.e., one-quarter of the period corresponding to each frequency after the rupture arrival. The rupture velocity changes on the fault plane and is 2.2 km/s from source 1 to source 2, and it changes to 2.6 km/s from source 2 to source 3. The amplitude has been normalized to the maximum value. The top, middle, and bottom panels show the image at the time when the wave has been radiated by source 1, source 2, and source 3, respectively. The stars represent source 1, source 2, and source 3, from left to right, respectively. (a), (b), (c) are results for the case = 2.2 km/s and = 2.6 km/s, = 1.5 km/s and = 2.6 km/s, and = 2.2 km/s and = 3.0 km/s, respectively. only signal waves but also noise correlated with the signal arrive at the array from different directions from the signal source, which reduces the accuracy of the incident direction estimation. Therefore, to suppress such signal correlation noise, it is desirable to configure the array with as many stations as possible, rather than the expected number of waves incident on the array +1. The array should be composed of as many stations as possible, with high correlation, matching polarity, and unification of whether the near-field term or far-field terms dominate the data.
In this study, in addition to the aforementioned conditions, we excluded stations with epicentral distances of >120 km because according to our preliminary study, the reflection wave from the Moho surface (SmS) reaches the stations with epicentral distances >120 km in the Kyushu region. Furthermore, to prevent the fault rupture image from being distorted by local site effects, only the underground stations of KiK-net (Aoi et al., 2004) were employed. Additionally, we visually checked the similarity of waveforms of adjacent stations to exclude records that were significantly different from other stations' waveforms, such as those strongly affected by site-specific effects. We also did not use stations near the causal faults (within approximately 30 km of the hypocenter), where there might be strong influences from the near-field term, intermediate term, and forward rupture directivity (Archuleta & Hartzell, 1981;Somerville et al., 1997). The use of stations at which the far-field term is dominant can prevent damage to the information on the rupture process, which is caused by the destructive interference among waveforms with poor correlations.

Strong Ground Motion Data
The acceleration waveforms from the underground KiK-net stations (Aoi et al., 2004) as shown in Figure 13 were corrected for the orientation of the seismometer installation and the offsets were removed. The waveforms were then integrated into the velocity using the methods reported by Kinoshita (1986). MUSIC (Schmidt, 1986) was used in the proposed method, which requires narrow-band waveforms. The, velocity waveforms were convoluted with Ricker wavelets using dominant frequencies of 0.25, 0.50, 1.00, 2.00, and 4.00 Hz to create waveforms for each frequency. We performed MUSIC and BP analyses with narrow-band waveforms obtained in this manner. Strong-motion records from KiK-net stations located within 120 km of the hypocenter were employed for the 2016 Kumamoto earthquake. Although KiK-net has both underground and surface stations, only the underground stations were used in this study to prevent the source image from being distorted by site amplification characteristics near the station. Based on visual confirmation of the observed waveforms and PS logging at the station, we selected the stations at which site effects were not significant. Even if the ground conditions of the station based on boring surveys were appropriate, if the observed records at the station were poorly correlated with the records of adjacent stations, the station was then excluded. Next, the selected stations that conformed to the above conditions were divided into three groups (FKOH, MYZH, and KGSH) according to the sign of the radiation coefficients of the SH waves. In each group, we visually verified the similarity of the waveforms, and stations with low similarity were excluded. Figure 12 shows the locations of the selected KiK-net stations. Figure 13 shows the observed north-south component velocity waveforms aligned according to the hypocentral distance, convoluted with the Ricker wavelet of the dominant frequency at 0.25, 0.50, 1.00, 2.00, and 4.00 Hz. If the P-wave portion was utilized in BP, only the waveforms from the arrival of the P-wave to the arrival of the S-wave was available, and we could not investigate the entire fault rupture process. Therefore, we used the S-wave portion; the first 0.5 s of the S-wave part was cosine-tapered while the P-wave portion before the S-wave arrival time was muted.

Fault Settings
Several fault models have been proposed for the 2016 Kumamoto earthquake based on surface displacements from ground surface surveys (e.g., Toda et al., 2016), aftershock distributions (e.g., Kato, Fukuda, et al., 2016;Kato, Nakamura, & Hiyama, 2016;Uchide et al., 2016;Yano & Matsubara, 2017), forward modeling of strong ground motion (e.g., Satoh (2017), and crustal deformation (e.g., Himematsu & Furuya, 2016). The common feature among these models is that the causal faults include both the Futagawa and Hinagu faults, with the Futagawa fault extending to Mt. Aso. In the BP analysis, specifying the fault plane in advance is not necessary; however, in this study, a fault model consisting of five fault planes was set based on Himematsu and Furuya (2016) and Yoshida et al. (2017) to improve the calculation efficiency. Ozawa et al. (2016) concluded that four fault planes are necessary to explain the crustal deformation that occurred during the 2016 Kumamoto earthquake based on InSAR and GNSS data. Based on the crustal deformation data from ALOS-2/PALSAR-2, Himematsu and Furuya (2016) identified the Hinagu fault, Futagawa fault, and an about 13-km-long fault approximately parallel to the Futagawa fault as the causal faults of the 2016 Kumamoto earthquake. They concluded that the slip on the fault plane, with a length of about 13 km, was a normal fault type, whereas the major faulting mechanism of the event was right-lateral strike-slip. The 2016 Kumamoto earthquake is the second example where such slip partitioning was observed (Toda et al., 2016 Figure 14 and Table 1 present the fault model assumed in this study.

Calculation of Travel Time
The BP method based on beamforming uses the travel time from the grid on the fault to the station while MUSIC uses the travel time difference between from the grid on the fault to the station and the from the reference point on the  fault to the station. The, the accuracy of the travel time data employed in BP or MUSIC significantly affects the accuracy of the results. In the BP of teleseismic records, the subsurface structure has a relatively minimal influence because the teleseismic waveforms propagate through the relatively homogeneous deep Earth. In contrast, in the BP using local strong-motion waveforms, seismic waves propagate in the shallow part of Earth characterized by strong heterogeneity. Therefore, realistic travel time data which reflects the heterogeneous velocity structure should be exploited to avoid distorting the image of the fault rupture due to the effects of wave propagation. Additionally, the causal faults of the 2016 Kumamoto earthquake were located in a volcanic area, where the underground structure is expected to be highly heterogeneous. Thus, in this study, the hypocenters of the 2016 Kumamoto earthquake and a series of seismic events that occurred around the causal faults of this earthquake were precisely determined using the 3-D velocity structure model reported in Matsubara et al. (2017), followed by the extraction of the observed travel times. By spatially interpolating the extracted 3-D observed travel times, we obtained the travel time from the computational grid on the fault plane of the 2016 Kumamoto earthquake to each station.
In the interpolation, data closer to the evaluation point should be given a larger weight. The weight was obtained based on the weighted average method shown in Equation 2, which has been used to generate equally spaced data via interpolations of heterogeneously distributed data (Inui et al., 2009;Shiono et al., 1985): where is the weight calculated according to the weighted average method from d, α, and S; d is the distance (km) between the data point and point where the interpolation is performed; and S and α are the parameters that determine how to weight the data according to the distance from the point where we aim to obtain the value.
Setting the values of S and α is necessary to perform the interpolation and minimize its error. Based on nth-fold cross-validation (Geisser, 1975), we chose the values of S and α to minimize the expected error. Generally, n should be at least 3 in the cross-validation; therefore, we conducted a 5th-fold cross-validation with n = 5. As a result, we found that the expected interpolation error became its minimum at S = 4.0 and α = 3.0; we therefore used these values to interpolate information on the relocated hypocenters and prepared the travel time data. The number of travel time data points for interpolation was 382562. However, if there were no more than three data points within a radius of 3 km of the grid for which we should obtain the travel time by interpolation, we projected the relocation data into the hypocentral distance versus the travel time space, interpolating the travel time in that space using Equation 2. In this case, we used the aforementioned values for S and α.

MUSIC Analysis
As described in Section 2.2, three arrays were configured from stations surrounding the source fault according to the sign of the radiation coefficient of the SH wave. As MUSIC uses relative travel times from a reference point, the absolute travel times from each grid point on the fault plane to each station were converted to relative travel times by subtracting the travel times from the hypocenter to each station. We employed the location of the hypocenter relocated with the 3-D velocity model reported in Matsubara et al. (2017) (130.764°E, 32.760°N, 12.5 km depth).
The MUSIC spectral image on the fault plane was calculated as described in Section 2.2. As seismic waveforms are non-stationary, the correlation matrix needs to be calculated adaptively; therefore, we estimated the correlation matrix sequentially and calculated the MUSIC spectrum. The waveforms used for MUSIC were normalized to the root mean square (RMS) amplitudes for the first 8.0 s of the S-wave component. Normalization by the average RMS amplitude should be more stable than that by the maximum amplitude, as the negative effects on the normalization owing to a sudden and abnormally large amplitude wave can be mitigated. This is a type of automatic gain control that compensates for the effects of amplification characteristics due to the propagation path and/or site. Both MUSIC and beamforming-based BP were performed with narrow-band waveforms, such that the amplification characteristics due to site effects and propagation path could be regarded as almost constant. Normalization by the RMS amplitude is expected to be effective as a correction for site and propagation effects.
The obtained MUSIC spectral images were binarized for each array according to the cumulative contribution ratio concept in PCA (see Section 2.3). Filtering was then performed based on the possible range of rupture velocity and maximum rise time. With reference to the velocity to set the timing of the first-time window trigger in the waveform inversion (2.4-3.0 km/s) by Asano and Iwata (2021), Kobayashi et al. (2017), Kubo et al. (2016), and Uchide et al. (2016), we set , , and as 2.0 km/s, 3.0 km/s, and 6.0 s, respectively, in Equation 1.

Beamforming-Based BP
BP was performed using all the arrays, as shown in Figure 12. To examine the entire rupture process, we executed a BP analysis using S-waves via a beamforming-based BP similar to that of , Kao and Shan (2004), and Krüger and Ohrnberger (2005). The Han window was not used in this study. BP was done using Equation 3 with an RMS envelope for the horizontal components: where ( ) is the radiation strength, is the location of the source point, is the time measured from rupture initiation, is the RMS amplitude of the S-wave, (t) and (t) are the north-south and east-west components, respectively, of the horizontal velocity record at the nth station, ( ) is the S-wave travel time from source point to the nth station, N is the number of stations. Here, represents the station correction, whose value acts to synchronize the initial parts of the S-wave of the waveforms from all stations when arranged according to the travel time from the hypocenter to each station. As can be seen from Equation 3, the radiation strength is a dimensionless quantity.
The amplitudes of the waveforms at each station were also normalized in BP. To calculate the RMS envelope, we synthesized the horizontal component with waveforms for both the north-south and east-west components. The sampling rate was 100 Hz. The optimal weights for each station were calculated according to Jakka et al. (2010) for each grid on the fault plane and used in BP. The BP results were multiplied by the binarized MUSIC spectral images and filtered using the information on the possible range of the rupture velocity and rise time.

Overall Features of the Rupture Process
Here, we summarize the overall features of the fault rupture process of the 2016 Kumamoto earthquake. Based on the distribution of the seismic radiation sources on the fault plane at each frequency (Figures 15-19), the overall image of the rupture of the 2016 Kumamoto earthquake shows that the rupture nucleated on the Hinagu fault (segment H in Figures 12, 14-19) and propagated in both the southwest and northeast directions. Significant seismic wave emissions on the Hinagu fault ended approximately 10.0 s after nucleation. The rupture on the Hinagu fault transferred to the Futagawa fault approximately 2-4 s after the earthquake. The major seismic wave radiation at depths from 0 to 15 km on the Futagawa fault occurred at approximately 6-14 s for 0.25-0.5 Hz and 6-18 s for 1-4 Hz after the origin time. This agrees with the spatiotemporal evolution of the fault rupture obtained by waveform inversion of the strong-motion records in Asano and Iwata (2021) and Kubo et al. (2016); these authors performed inversion in the frequency band of 0.05-0.50 Hz and 0.05-1.0 Hz, respectively. In this study, BP at frequencies of 0.25, 0.5, 1, 2, and 4 Hz was done, hence an approximate comparison of the spatiotemporal evolution of the fault rupture is possible. We also found that the F4 segment parallel to the Futagawa fault, estimated from crustal displacement data in Himematsu and Furuya (2016), ruptured at the same time as the onset of the major slip on the Futagawa fault (approximately 6-12 s after the earthquake origin time) (Figures 15-19). Additionally, the spatiotemporal distribution of the seismic radiation sources varied depending on the frequency.

Rupture Process at Each Frequency
We analyzed the distribution of seismic wave radiation sources at 0.25, 0.5, 1, 2, and 4 Hz on the fault planes. In view of the results at 0.25 Hz in Figure 15, the seismic wave emission source evolved bilaterally from the hypocenter to the southwest and northeast on the Hinagu fault. The rupture transferred from the Hinagu fault to  Table 1 and Figures 12 and 14. (b) The raw seismic radiation by beamforming. Neither a multiplication of MUSIC spectral image nor rupture velocity filter has been applied.
the Futagawa fault 2-4 s after the earthquake initiation as mentioned above. This is consistent with the spatiotemporal evolution of fault rupture obtained by waveform inversion in Kubo et al. (2016) (Figure 21). In our study, a significant seismic wave emission is seen on the Futagawa fault at depths <15 km 6-12 s after the earthquake occurrence. This also agrees with the inversion results of Asano and Iwata (2021) and Kubo et al. (2016). A substantial seismic wave radiation was also detected at a depth of 5-20 km on the Hinagu fault 4-10 s after earthquake onset. In contrast, seismic wave emission at 0.5 Hz ( Figure 16) and 1 Hz (Figure 17) in this area is not significant. This illustrates that the distribution of the seismic radiation area on the fault plane significantly differs depending on its frequency.
In Figure 16, at a frequency of 0.5 Hz, seismic wave emission can be observed across a wide area at depths of 0-15 km on the Futagawa fault 6-12 s after the origin time. Seismic wave emission is also observed at depths of 5-15 km on the Hinagu fault until approximately 10.0 s after the earthquake onset. The range of seismic wave emission shows that the depth distribution of the seismic wave emission source at a frequency of 0.5 Hz covers the area at depths from 0 to 20 km on the fault plane; however, most of the sources are located in areas shallower than 15 km. The major seismic radiation at 1 Hz ( Figure 17) also expands from depths of 0-20 km; these areas radiated significant seismic waves 4-18 s after the earthquake initiation. The seismic wave emission sources at 1 Hz are relatively more widely distributed in space compared with those at 0.25 and 0.5 Hz. This also indicates a difference in the distribution of the seismic radiation sources depending on the frequency, which has been reported in previous studies for some earthquakes (e.g., Nakahara, 2008).
At frequencies of 2 ( Figure 18) and 4 Hz (Figure 19), the seismic radiation areas are distributed from depths of 0-20 km; these areas radiated seismic waves at approximately 4-18 s after the origin time. The area of each significant seismic radiation source is small compared with the seismic radiation areas at 0.25, 0.5, and 1 Hz. The seismic radiation source distribution at 2 and 4 Hz is wider than those at 0.25, 0.5, and 1 Hz. The higher the frequency, the smaller the size of each seismic radiation source and the wider the spatial distribution. Such a frequency-dependent feature in the seismic radiation distribution has implications for the establishment of the seismic wave radiation sources on the fault plane for strong broadband ground motion predictions. Figure 20 shows the normalized accumulated seismic radiation on the fault obtained for each frequency. The seismic radiation sources are mainly distributed on the Futagawa fault from depths of 0-20 km; most of the accumulated seismic radiation sources expand at depths from 0 to 15 km. They are contained in the large fault slip area obtained by Asano and Iwata (2021) and Kubo et al. (2016) (Figure 21) with the waveform inversion of strong-motion records. This study targeted different frequency bands from those studies, making direct and rigorous comparisons difficult. However, a rough comparison could be made for the frequency band from 0.25 to 1 Hz. At 0.5, 1, 2, and 4 Hz, large seismic radiation was observed near the ground surface of the Futagawa fault at a horizontal distance of 0-25 km from the hypocenter. In this area, ground surface surveys (e.g., Toda et al., 2016) and SAR image analyses (e.g., Himematsu & Furuya, 2016) have reported large displacements; our results showed that large seismic radiation occurred in these large slip areas. However, no significant seismic radiation was detected at 0.25 Hz near the ground surface in this area. The BP method used in this study is different from that used in Festa and Zollo (2006). They directly treated displacement on the fault as an unknown parameter while determining the displacement by linear inversion of the high-frequency contents in the observed records, with theoretical Green's functions. In contrast, our BP method is based on beamforming where waveforms recorded at stations are stacked by adding appropriate time delays. We used only velocity waveforms of stations distant from faults, where the far-field term is dominant. Figure 22 is the histograms of the number of stations used in this study and the shortest distance between the station and each fault. We can see that all the stations are distant from faults and waveforms are expected that far-field term is dominant. Static displacement caused by earthquakes is contained in near-and intermediate-field terms (Aki & Richards, 2002) in the velocity waveforms. Since we used only the waveforms in which the far-field term is dominant, little static displacement is supposed to be included in the waveforms. There is a possibility that this caused no detection of static fault displacement in our results even at a low frequency of 0.25 Hz. Figure 23 illustrates the slip distribution on the fault plane obtained by Himematsu and Furuya (2016) by inverting the ground surface displacement obtained from the pixel-offset analysis of SAR images (note that the depth of the upper edges of all fault segments is set to zero). This figure shows that there is large fault displacement in a relatively shallow northeast part of the Futagawa fault, where no seismic emission sources were detected in our results at 0.25 Hz, similar to the results of Asano and Iwata (2021) and Kubo et al. (2016). As mentioned above, static displacement is included in the near-and intermediate-field terms of velocity waveforms and is almost not included in the waveforms we used for BP in which the far-field term is dominant. Thus, our results for shallow part at 0.25 Hz are not consistent with those by Asano and Iwata (2021), Kubo et al. (2016), and Himematsu and Furuya (2016). Asano and Iwata (2021) performed an inversion of the 2016 Kumamoto earthquake using waveform records from multiple stations very close to the fault where the near-field term predominates. They quantitatively investigated where on the fault plane the prominent pulses with a period of 0.25-0.3 Hz were observed, and they found that a significant portion of them originated at depths shallower than 3 km to shallow depths of the fault, in harmony with conclusions by Kobayashi et al. (2017). This means, conversely, that this displacement at the shallow part of the fault, which is the main cause of seismic radiation at around 0.3-0.25 Hz observed only in the area very close to the fault, cannot be detected from waveforms consisting only of the far-field term.

Seismic Wave Radiation Distribution and Waveform and SAR Inversion Results
As mentioned before, we used only velocity waveforms in which the far-field term is prevailing. The velocity waveforms of the far-field term are proportional to the slip acceleration of the fault (e.g., Aki & Richards, 2002). Therefore, as also stated by Yao et al. (2013), looking at the seismic radiation intensity from a fault by stacks of velocity waveforms consisting of far-field terms correspond to the slip acceleration on the fault. In general, since the waveform inversion analyzes the slip or slip velocity of the fault, the BP captures seismic radiation with a higher frequency component emphasized than the waveform or crustal deformation data inversion. Even if BP and waveform inversion uses the same bandwidth of data, they show different phenomena, and the seismic source distribution on the fault plane detected by BP does not necessarily match the displacement of the fault itself as Figure 20. Normalized accumulated seismic radiation on the fault at each frequency. The radiation strengths were normalized to the maximum value on the fault plane at 0.25 Hz. The approximate maximum radiation strength normalized by that at 0.25 Hz is shown on the right side of each panel. Open circles depict the hypocenters of the aftershocks that occurred from 04/2016 to 12/2016 with a distance from the fault plane within 5 km. obtained through waveform inversion or tectonic data inversion. Yao et al. (2013) also concluded that BP results cannot be directly related to fault slips. Additionally, Walker and Shearer (2009) reported that the seismic wave radiation source distribution obtained by BP of teleseismic P-wave velocity waveforms correlates poorly with the seismic moment distribution obtained by waveform inversion, and relatively better with the seismic moment rate. In accordance with their results, the area detected by our BP method corresponds to an area with a large slip rate, rather than an area with large displacement and small slip velocity. Compared with the slip velocity in the region with a large accumulated seismic radiation in this study, the estimated slip velocities in Asano and Iwata (2021) and Kubo et al. (2016) was also large. The results of this study agree with those of Walker and Shearer (2009). At a horizontal distance of approximately 0-10 km and 25-35 km from the hypocenter, there were only moderate seismic wave emissions near the ground surface at all frequencies. This roughly corresponds to the location where the fault slip velocity is low with a moderate slip amount, as obtained via waveform inversion in Asano and Iwata (2021) (Figure 21). Himematsu and Furuya (2016) reported that the F4 segment (a small fault parallel to the Futagawa fault; the rightmost segment in each panel) had a large slip that reached the ground surface from an inversion analysis of pixel-offset data of SAR images (Tobita et al., 2001). The accumulated seismic radiation from 0.5 to 4 Hz, as shown in Figure 19, is significant, which may be due to the detection of a large slip acceleration on this small fault during the earthquake. In our study, high-frequency seismic radiation at 2 and 4 Hz was found in shallow parts of faults F2, F3, and F4. The possible cause of this may be that the BP tends to catch seismic radiation images with an exaggeration of the higher frequency component, as discussed here.

Accumulated Seismic Radiation Distribution at Each Frequency
Although the accumulated seismic wave radiation distributions at 0.25, 0.5, and 1 Hz slightly overlap each other, the largest seismic radiation areas at each frequency are different. On the Futagawa fault, at depths deeper than 5 km in the area approximately 30-50 km from the hypocenter, both slip and slip velocity are large in Kubo et al. (2016) and Asano and Iwata (2021) (Figure 21). In their results, although the slip amount is large, the slip velocity is small in the area shallower than 5 km in depth. The static displacement in this region was not detected in this study because we used only waveforms dominated by the far-field term. Excluding this shallow area, the accumulated seismic wave radiation at 0.25 Hz is also significant in the area where both slip and slip velocity inferred by Kubo et al. (2016) and Asano and Iwata (2021) are large. Substantial seismic radiation areas at 1-4 Hz are distributed around those at 0.25 Hz. In this area, the higher frequency (1-4 Hz) seismic radiation sources are almost at the rim of the lower frequency (0.25 Hz) seismic radiation sources, and their distribution is complementary to each other. Similar results were reported by Nakahara (2008). He also found the fact that the distributions of the low-and high-frequency seismic radiation sources are complementary in some earthquakes (e.g., Sanriku-oki earthquake (M7.6) in 1994, Northern Iwate earthquake (M5.9) in 1998, Chi-Chi Taiwan earthquake (M7.6) in 1999, etc.) which is in agreement with the observations of the present study. Such high-frequency radiation at the periphery of the large moment release area estimated by waveform inversion is also reported in  for the 1993 Kushiro-Oki earthquake (M7.8) and  for the 1995 Hyogo-ken Nanbu earthquake by the envelope inversion of acceleration seismograms. The existence of the high-frequency radiation around the low-frequency radiation region may be attributed to (a) the concertation of high-frequency seismic radiation due to the slip rate and stress at the fault edge (Madariaga, 1977) and (b) the emission of high-frequency seismic radiation that occurs with the arrest of slip (associated with the stopping phase) (Madariaga, 1976). In the slip rate functions reported in Asano and Iwata (2021) and Kubo et al. (2016) (Figure 21), the duration of the slip rate function at a high amplitude is longer at the center and shorter at the rim of the large slip area, characterized by a shape similar to that of the slip rate function analytically deduced by Kostrov (1964). This indicates that lower (higher) frequency seismic waves radiated at the center (rim) of the large slip area. The distribution of seismic radiation obtained in our study is consistent with the conclusions about the mechanism of the high-frequency radiation by Madariaga (1976Madariaga ( , 1977 and the slip rate function put forth by Asano and Iwata (2021) and Kubo et al. (2016). Our results show that the distribution of the seismic wave radiation sources depends on the frequency. It might suggest that the placement of a strong seismic radiation area should have the frequency dependence for seismic motion prediction.

Accumulated Seismic Radiation Distribution and Strong-Motion Generation Area (SMGA) Model
We compared the distribution of the accumulated seismic radiation in Figure 20 with the fault model developed by Satoh (2017, as shown in Figure 24. Satoh (2017 identified five SMGAs, SMGA1-5, to reproduce the 2016 Kumamoto earthquake ground motion by forward modeling using the empirical Green's functions method. She modeled the SMGAs as areas with a large slip and stress drop that radiated strong motion. Figure 24 a shows the seismic source model proposed by Satoh (2017; as the source model was developed to reproduce seismic motion from 0.2 to 10 Hz, it might not be directly comparable to the seismic wave radiation at each frequency obtained in this study; however, a rough comparison could be possible. Comparing Figures 20 (or 24b) and 24a, it is evident that the seismic radiation source areas detected in this study generally agree with the location of the SMGAs identified by Satoh (2017); SMGA1 (at 0.25 Hz), SMGA2 (at 0.25, 0.5, and 1 Hz), SMGA3 (at 1, 2, and 4 Hz), and SMGA4 (at 0.5, 1, 2, and 4 Hz) were detected as strong seismic radiation areas. For SMGA5, only moderate seismic radiation is detected (at 1 Hz), which likely originated from the small slip rate in this area (the top left part of the slip rate functions on the fault as shown in Figures 21a and 21b). Overall, the location of the SMGAs in Satoh (2017 and the large accumulated seismic radiation areas obtained in this study naturally coincide. It shows that the method proposed in this study effectively detects SMGAs.

Distribution of Seismic Wave Radiation and Aftershocks
For many earthquakes, the distribution of the fault slip and aftershocks are complementary owing to the stress released by the large slip of the mainshock (e.g., Doser & Kanamori, 1986;Hartzell & Heaton, 1986;Mendoza & Hartzell, 1988). Based on Figure 20, some aftershocks including a moderate event occurred in significant seismic radiation areas at 0.25 and 0.5 Hz; it is not clear whether the aftershock distribution is complementary to the seismic radiation distribution. Although the aftershock activity was highly dynamic adjacent to the causal fault    Uchide et al. (2016) and Yano and Matsubara (2017) reported the presence of an inactive aftershock area, which is clearly visible in Figure 20 (the distance along the strike is approximately 35-45 km and deeper than about 10 km). These authors discuss the reason(s) for the absence of seismic activity in this region. They indicate that slip during the main shock may have released stress and reduced the number of aftershocks. Kato, Nakamura, and Hiyama (2016) also pointed out that fault slip or thermal structures may have produced areas of low aftershock activity. The source fault of the 2016 Kumamoto earthquake passed through the interior of Mt. Aso, an active volcano. Geophysical investigations such as seismic tomography (e.g., Sudo & Kong, 2001;Wang et al., 2018), gravity surveys , and seismic attenuation tomography (Komatsu et al., 2017), have reported the heterogeneity of underground structure in and around Mt. Aso. Their results indicate that this region is heterogeneous in terms of its dynamics and physical properties. The influence of these heterogeneities on the source process of the 2016 Kumamoto earthquake has also been discussed, including the possibility that a magma pool terminated the fault rupture (e.g., Miyakawa et al., 2016;Yagi et al., 2016), or that there was a material barrier due to the presence of the Aso caldera (Yue et al., 2017). As shown in Figure 20, the BP analysis conducted in this study did not indicate any significant seismic emission at distance of 35-40 km along the strike and at depths of more than about 10 km. In this area, the slip and slip velocity are also small as reported by Asano and Iwata (2021) (distance along strike around −25 km in Figure 21a) and Kubo et al. (2016) (the lower part on the left side of the fault in Figure 21b with the smallest slip velocity). This may be due to the high temperature and non-brittle state of this region, which renders the emission of significant seismic waves difficult and possibly also suppresses the aftershock activity. In contrast, the region along the strike around a distance of 40-45 km and a depth of more than about 10 km shows moderate seismic radiation as shown in Figure 20. The slip velocity function obtained by Kubo et al. (2016), as shown in Figure 21b, also has a moderate value around this area (the leftmost part of the fault). Considering that the BP method detects the rate of slip rather than the amount of slip itself, the magnitudes of slip and slip velocity in this region may have been moderate during the main shock, which in turn may have led to the lack of aftershocks. However, only very small slip values and slip velocities exist in the results reported by Asano and Iwata (2021) as shown in Figure 21b. Hence, it remains inconclusive whether the suppression of the aftershock activity in this region was due to the release of stress caused by the mainshock or due to the subsurface structure and its physical properties. As for the area along the strike at a distance of 35-40 km, and a depth greater than ∼10 km is likely to be non-brittle and characterized by high temperatures, given the lack of aftershocks and the fact that little seismic wave radiation or fault slip occurred during the earthquake.

Summary of Results and Discussion
Based on the heterogeneity in the seismic radiation and aftershock distributions (Figure 20), the 2016 Kumamoto earthquake occurred in a highly heterogeneous region in terms of its mechanical and physical properties. The 2008 Iwate-Miyagi Nairiku earthquake, which occurred along an active fault in a volcanic zone, was controlled by thermally and physically heterogeneous structures from its onset throughout the fault rupture process (Ando & Okuyama, 2010;Oshima & Takenaka, 2022). The mechanical and physical heterogeneities peculiar to volcanic regions are important components in strong-motion predictions and hazard assessments for earthquakes resulting from active faults in these regions. We showed that a detailed rupture process can be revealed even in volcanic areas by applying the proposed method with accurate 3-D travel time data. Using this method, the rupture process of earthquakes in volcanic regions can be elucidated and the data accumulated can be used to predict strong earthquake motion in volcanic regions. Additionally, the present method can detect frequency-dependent seismic radiation source distributions. The fault models used for broadband seismic motion prediction should consider the differences in the seismic radiation source distribution by frequency, as pointed out in previous studies, as well as in this study. Particularly, the detailed distribution of seismic radiation sources on the fault plane must be accurately modeled for strong ground motion predictions near faults, in which the near-field and intermediate terms are dominant. The proposed method can be used to investigate the fault rupture process at different frequencies with a high resolution while suppressing artificial noise; hence, it may be useful in understanding and modeling the details of fault rupture for predictions of earthquake ground motion adjacent to a fault.

Conclusions
By combining BP based on beamforming and MUSIC, we developed a BP method for imaging fault ruptures with noise suppression and high resolution. The stations surrounding the fault were grouped into arrays according to the sign of the radiation coefficients of the SH wave. MUSIC spectral images were then obtained from each array and were binarized according to the accumulated contribution ratio in PCA. The MUSIC spectral images from the arrays were multiplied, and the seismic radiation areas on the fault planes were extracted after filtering based on the possible range of rupture velocity and assumed maximum rise time. The extracted radiation regions were given a value of 1 while the remaining regions on the fault plane were set as 0 for use as weights. These weights were multiplied by the BP image, which was obtained using all the arrays to capture a high-resolution rupture image. By performing this process at each time step, high-resolution BP analysis with noise suppression was possible. To improve the accuracy of MUSIC and beamforming, we interpolated the travel time data from a relocated event catalog using the 3-D seismic wave velocity structure model by Matsubara et al. (2017). The optimal values for the control parameters in the interpolation were determined via 5-th fold cross-validation.
We tested the present method on a rupture with varying the velocity on a fault plane, using sine waves with noise, theoretical waveforms with radiation coefficients computed using a 1D velocity structure model. In the tests, we also examined the effect of the assumed rupture velocity on the results. The results show that noncorrelated random noise performs well for rupture velocities varying on the fault plane up to a signal-to-noise ratio of −16 dB. Tests considering the radiation coefficient and 1D subsurface model confirmed the validity of the proposed method, since the large radiation strength at seismic sources was detected, although the performance of beamforming and MUSIC was degraded by the later phases. It was also shown that the assumed rupture velocity controls the results of the proposed method and that it is necessary to set the value as close as possible to the average rupture velocity to minimize artifacts as much as possible. It is also found that the proposed method is robust concerning the assumed rupture velocity as long as it does not deviate too much from the average rupture velocity of the target event.
We applied the present method to the strong motion records of the 2016 Kumamoto earthquake. Images of the fault rupture at 0.25, 0.5, 1, 2, and 4 Hz were captured, showing that the rupture initially propagated bilaterally from the hypocenter. The rupture then transferred to the Hinagu fault and ended 10.0 s after the nucleation. Finally, the rupture propagated to the Futagawa fault 2-4 s after the earthquake onset. Major seismic wave radiation, at depths of 0-15 km on the Futagawa fault, occurred 6-14 s for 0.25-0.5 Hz and 6-18 s for 1-4 Hz after the origin time. These results were consistent with the spatiotemporal evolution of the fault rupture estimated via the waveform inversion of strong-motion records by Asano and Iwata (2021) and Kubo et al. (2016). The normal fault-type slip on the F4 segment parallel to the Futagawa fault, as inferred from the crustal deformation data by Himematsu and Furuya (2016), occurred at the same time as the main slip on the Futagawa fault (6-12 s after the origin time).
The distribution of the accumulated seismic radiation obtained in this study correlates to the fault slip velocity than the fault slip amount. The area where accumulated seismic radiation was large tended to coincide with the high slip rate area obtained from waveform inversion results. This may be because only the velocity waveform, in which the far-field term predominates, was stacked in our BP process, resulting in imaging of the slip acceleration of the fault. It means that we are looking at a source image in which relatively higher frequency components are emphasized than in the waveform inversion. This may be the possible cause of the seismic radiation at 2.0 and 4.0 Hz detected in the shallow part of the Futagawa fault in this study. In the shallow part of the Futagawa fault, where a large fault slip had been identified from waveform inversion and geodetic data inversion studies. However, no cumulative seismic radiation at 0.25 Hz was observed in this study. This is because static displacement is generated by the near-and intermediate-field term part of the waveform, while only the velocity waveform dominated by the far-field term was used in this study, and the near-and intermediate-field terms mainly caused by movement of the shallow portion of fault plane was not considered in out BP.
Our BP images detected the area of large slip and slip rate around 30-50 km from the hypocenter on the Futagawa fault, where the spatiotemporal distribution of the seismic radiation sources significantly varied according to the frequency; the spatial distributions of the low-and high-frequency seismic sources were complementary. The high-frequency (1-4 Hz) seismic radiation sources occurred at the rim of the low-frequency (0.25 Hz) ones as certain previous earthquakes reported in previous studies. This could be naturally understood in the context of the high-frequency seismic radiation mechanism: the high slip rate and stress concentration at the fault edge and the high-frequency seismic radiation at the stopping phase associated with slip arrest. It suggests that future ground motion predictions should take into consideration the setting of strong seismic radiation areas and their frequency dependence. The distribution of the SMGAs modeled by empirical Green's function method using the strong motion data was also consistent with the distribution of the large cumulative seismic radiation area obtained in this study. This reinforces the aforementioned conclusion that the large cumulative seismic radiation area identified in this study corresponds to regions of large fault slip rates, as SMGAs are regions that generate strong ground motion and are characterized by large fault slip rates.
The proposed method effectively suppressed artifacts in the BP images, which is robust to noise, and give high-resolution results. This method can provide a reference to understand and model the details of seismic radiation sources, which is important for strong ground motion prediction in near-fault areas.

Data Availability Statement
We employed waveform data from KiK-net (Aoi et al., 2004) provided by NIED. The strong ground motion data is available at KiK-net's website: www.kik.bosai.go.jp (last accessed October 2022). We used the F-net solution (Fukuyama, 1996;Fukuyama et al., 1998) from NIED for the moment tensor solution of the 2016 Kumamoto earthquake. Generic Mapping Tools (Wessel et al., 2013) were used to create several figures in this paper.