Estimating Rupture Front of Large Earthquakes Using a Novel Multi-Array Back-Projection Method

A ruptured front obtained from high-frequency energy radiation is the key to understand the complex source. It is commonly observed that rupture fronts derived from different arrays often show some variations due to the obvious difference of the positioning accuracy of the far-field array between the azimuth and the epicentral distance. We developed a new multi-array back-projection method based on the classical back-projection method and applied the method to the 2015 M W7.8 Nepal earthquake. The back azimuth information with small error is separated from the classical back-projection results, and the azimuth intersection of multiple arrays is used to obtain more accurate spatial and temporal distribution information of the source rupture fronts.


INTRODUCTION
The array back-projection (BP) method was first successfully used to image the source rupture process of the 2004 Sumatra earthquake (Krüger and Ohrnberger, 2005a;. In this method, the array signal processing technology is used to analyze the seismic waves recorded at the seismic network, so as to image the position of the source rupture front and obtain the spatial-temporal distribution image of the source rupture process. The back-projection method has been successfully applied to several major global earthquakes Kiser and Ishii, 2011;Meng et al., 2011;Wang and Mori, 2011;Zhang et al., 2011;Koper et al., 2012;Yao et al., 2012;Fan and Shearer, 2015;Meng et al., 2016;Wang and Mori, 2016;Liu et al., 2017). At present, it has become a widely used method for source process imaging. This method differs from other approaches because it requires less of a hypothesis of the source model, so it can estimate the rupture speed independently (Wang et al., 2016a). It uses high-frequency signals, so it can give more detail about the source rupture process. Moreover, this method substantially reduces the effort involved in program design and computing time.
In essence, the back-projection method is a continuous relative positioning technology. As the seismic array is located on the side of the seismic source, the azimuth coverage is smaller than that of the seismic network, and there are systematic errors in positioning, which are reflected in the systematic deviation of the position of the rupture point and smeared in the space domain (Krüger and Ohrnberger, 2005b). In addition, the positioning error due to the three-dimensional earth structure makes the results between different arrays inconsistent. Therefore, the initial rupture position information of the main shock is used to correct for the travel time anomalies caused by the three-position velocity structure.
With the propagating of the rupture front, the hypocenter correction gradually becomes ineffective for earthquakes with a large extent. It is possible to consider treating the global network as a platform array (Walker et al., 2005;Xu et al., 2009) to increase the azimuth coverage and reduce the systematic error of array positioning. However, due to the directional nature of seismic energy radiation, the correlation between the highfrequency signals recorded by stations in different orientations is weakened, leading to a sharp decrease in the number of available stations.
Another option is introducing the aftershock position information to correct the back-projection imaging position in the middle and late stages of the rupture process (Ishii et al., 2007;Palo et al., 2014;Meng et al., 2016;Qin and Yao, 2017). Wang et al. (2016b) proposed the image deconvolution back-projection method. In Wang's method, moderate aftershock is introduced as the reference event, and the array response can be removed by deconvolution to obtain a sharp image of the high-frequency energy radiation. However, it cannot eliminate the swimming artifact in the late period of the rupture process. The slowness correction method proposed by Meng et al. (2016) can better solve the problem of inconsistent results between seismic arrays, and a small number of aftershocks must be used.
The common limitation of location correction methods that operate by means of aftershocks is that they strongly depend on the accuracy of the location results of aftershocks and the happening of large aftershocks. Many studies suggest that the uncertainties in aftershock locations are much larger than the uncertainty of main shocks, and it is difficult to obtain timely and reliable locations of aftershocks. Moreover, the distribution of aftershocks is not necessarily uniform, and the location of large aftershocks may not be on the rupture trace of the main shock, so the imaging results of multiple arrays cannot be given in a timely and accurate manner after the earthquake.
In order to effectively utilize the azimuth coverage capability of multiple arrays, Xu et al. (2009) applied teleseismic multiple array back-projection (Huang et al., 2012;Kiser and Ishii, 2012;Zhang et al., 2016;Qin and Yao, 2017) to study the rupture process. They back-projected the P wave to the source area, respectively, for each array and then superimposed the energy in the same time window for finding the position of the maximum superposition energy in each time window, so as to obtain the spatiotemporal process of the earthquake source rupture. However, due to the radiation pattern and the heterogeneities in the ray paths, it leads to different locations of each array back-projection. In the process of multi-array stacking, the array with the largest beam energy value usually plays a major role. Although the weight factor can be introduced, the positioning information of multiple arrays cannot be integrated effectively.
According to the published results, there are still great differences in the source kinematics parameters of the same earthquake obtained using different seismic arrays. Based on the array response function and the existing research results, we analyze the positioning errors of the array back-projection method in the radial and tangential orthogonal directions and find that the locating errors in the two directions have great differences and dynamic changes. We use this property to develop a multi-array cross back-projection method (MCBP) based on the classical array back-projection method and test the effect of this method by taking the 2015 Nepal M W 7.8 earthquake as an example. The method separates the back azimuth information with small error from the results of the classical array backprojection method and crosses the back azimuth to obtain the spatial and temporal distribution of the high-frequency radiation sources.

Location Bias in Classical Back-Projection of the 2015 Nepal-Gorkha Earthquake
The 2015 Nepal-Gorkha M W 7.8 earthquake caused more than 9,000 casualties, making it the most devastating and disastrous earthquake to hit Nepal since the 1934 Nepal-Bihar earthquake. The abundant array records of the Nepal-Gorkha earthquake and sufficient research provided an ideal case for understanding and testing the MCBP method.
We downloaded 649 broadband vertical component seismograms with an epicentral distance ranging from 30 to 90°from the IRIS DMC. By calculating the correlation coefficient 15 s after the initial P wave, 292 seismic records with a correlation coefficient greater than 0.7 were selected, constructing five regional arrays ( Figure 1). They are the Europe array (EU, 65 stations), the Alaska array (AL, 90 stations), the Japan array (JP, 28 stations, most of them are in Japan), the Australia array (AU, 90 stations), and the southern Africa array (AF, 19 stations). They are roughly five-pointed stars.
The focal area (27°N ∼ 29°N, 84°E ∼ 87°E) was divided into a 0.05°× 0.05°grid with a sliding time window of 15 s and a step length of 1 s. Considering that the back-projection method can suppress the noise from outside the seismic source, we use the back-projection method to choose the maximum available frequency. To make full use of the high-frequency signal, the high-frequency information of 0.3-10 Hz was used. The location of the epicenter was determined using USGS (28.2305°N, 84.7314°E). Based on the AK135 velocity model, the relative back-projection method (Du et al., 2009;Zhang and Ge, 2010) was used to back-project the waveform onto the horizontal fault surface with a depth of 8 km.
The high-frequency radiation sources ( Figure 2) obtained by the five-array imaging have good consistency in the first 40 s, all extending southeast from the initial rupture point, but there are significant differences between each array.
In order to analyze the positioning differences between different arrays, we calculated the extension response functions of the generalized array based on spherical waves (modified from Xu et al., 2009).
where θ is the longitude, ϕ is the latitude, h is the depth, f is the frequency, N is the number of seismometers, t(θ, ϕ, h) j is the travel time from point (θ, ϕ, h) to the jth seismometer, and t(hypocenter) j is the travel time from the hypocenter to the jth seismometer.
The extension array response functions of the five arrays are shown in Figure 3, and we found that most array response functions are ellipses. According to the results of the back-projection method published so far, the distribution of back-projection stacking energy is generally elliptical. This distribution characteristic indicates that the error of the position of the high-frequency radiation sources obtained by the back-projection of the array is large in a certain direction, which may lead to a sudden jump in the position of the adjacent radiation sources. Krüger back-projected the GRSN array and the FNET array, respectively, to obtain the rupture process of the Sumatra earthquake (Krüger and Ohrnberger, 2005b). The back-projection results of the GRSN array showed a jump in the second half of the rupture, while the results of the FNET array showed an obvious jump in the first half of the rupture.
Array response function is a good indicator of the array resolution, but it ignores the heterogeneities in the ray paths, S/N ratios of the waveforms, and instrument qualities. To show a comparison for the results derived from different arrays, we can back-project several aftershocks using the travel time corrections derived from the main shock. There are researchers who have done this (Meng et al., 2016;Wang and Mori, 2016). In both Wang and Meng's respective studies, the aftershocks were backprojected. They all found that the aftershock locations inferred  can be thought to be consistent. To synthesize the above analysis, we think that the resolution of back-projection is not isotropic, and the spatial bias varies with the position of the rupture front. The long axis of the error ellipse is always toward the array. The two orthogonal directions have different physical significances, the radial corresponding to the epicentral distance of the source and the tangential corresponding to the back azimuth. The resolution analysis shows that the error of the back-projection method is larger in measuring the epicentral distance and smaller in measuring the back azimuth. Based on this characteristic, we design a multi-array back-projection method.

Multi-Array Cross Back-Projection Method
We decompose the results from the classical array backprojection method into two orthogonal directions, radial and tangential. The error analysis above shows that the error of the seismic array in estimating the back azimuth is less than that in estimating the epicentral distance. We can use the back azimuth alone. Considering that only one inverse azimuth can be determined by a single array, theoretically, at least two arrays of different azimuths are needed to determine the location of the focal point by the back azimuth component. It can be divided into several sub-arrays for the global network or the regional network, which is smaller than the aperture of the whole array, so as to avoid the influence of the source radiation pattern and the Doppler effect of rupture propagation to the greatest extent.
In order to separate the back azimuth from the results of the classical back-projection method, we connect from the reference station at the center of the array to the highfrequency radiation sources determined by the backprojection method (Figure 4). The location error of the epicenter distance is large, so we can assume that the true location of the subevent within each time window is not determined, but it is constrained on this ray. If there are two or more arrays located in different directions of the epicenter, multiple rays can be determined. The intersection point of the rays can determine the best estimate of the locations of the high-frequency radiation sources.  This method needs to ensure that the signals in the same time window of different arrays are homologous. The source of a large earthquake is a moving source with unilateral or bilateral expansion, and the magnitude of its expansion velocity is the same order of magnitude as the magnitude of the P wave propagation velocity, which must lead to the observation in different azimuth arrays, and the duration of the rupture is significantly different. In front of the rupture, the wave train compresses, and the duration becomes shorter. Conversely, in the direction away from the rupture, the wave train is elongated. However, the back-projection generally uses uniform sliding time windows, and the waveform in the same sliding window on different arrays may come from different sub-sources. In order to ensure the signal homology, we need to modify the backprojection results of each array.
It is assumed that the travel time of the P wave radiated by the epicenter to the array is τ 0 . After the rupture process lasts for t r , the travel time of radiated seismic waves from a rupture point to the array is τ 1 . Then the time difference between the arrival of the subsequent wave train and the first arrival is t r ′ , as calculated below: and the derivation gives us the following: This formula indicates that we can deduct the time difference caused by the change of the position of the rupture point from the apparent rupture time to obtain the rupture duration at the source. The time difference t r ′ can be obtained using the relative back-projection method (Zhang and Ge, 2010). The travel time difference (τ 1 − τ 0 ) can be calculated according to the AK135 velocity model. Since the time interval of the modified rupture points is no longer uniform, it is necessary to obtain the rupture process with equal time intervals through interpolation processing.
Ideally, the point of intersection is a geometric point. If there are more than two arrays, it is difficult to get an ideal geometric point directly. We can fit the equation of the ray and establish the linear equations. Then we use the least square method to solve the linear equations and calculate the position of each intersection point.

Numerical Test
In order to evaluate the feasibility of the new method and its ability to resolve the spatiotemporal position of the rupture point, we performed a numerical test using the theoretical seismogram calculated using QSSP . We take the epicenter of the Nepal-Gorkha earthquake (28.23°N, 84.73°E) as the initial rupture point. It is assumed that 10 sub-sources with a thrust mechanism (strike 287°, dip 6°, slip 96°) rupture occur in succession at a velocity of 2.5 km/s. We set an inflection point in the middle of the rupture, and the focal depth increased from 8.2 to 18.2 km. The rising time of the source function is set as 1, 1.5, 3. 5, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, and 1 s, respectively. The scalar seismic moment of each sub-source is proportional to its rising time, and the sum of the seismic moments is equal to the scalar seismic moment of the Nepal earthquake (8.39 × 1020N·m). The total rupture length is 200 km, the duration is 81 s, and the seismic wave sampling rate is 5 Hz. The AK135 reference model was used in the calculation. There are five arrays in the global network that record analog data.
First, we use the standard single array technique to calculate the position of the rupture sources ( Figure 5). We can see that the imaging results of each array are significantly different, and most arrays fail to accurately obtain the true input source location. Next, based on the imaging results, we applied our new method to obtain the imaging results of multiple arrays. The results ( Figure 6) show that the spatiotemporal location of the rupture point obtained using the new method is basically the same as the real location. Due to the smoothing effect of the sliding time window, the new method can yield a series of rupture points. These rupture points are basically located on the curve formed by the 10 sub-sources. Especially at the inflection point, the direction change of rupture propagation calculated using the new method is in good agreement with the actual direction change. It is fully proved that the new method has high resolution for earthquakes where the direction of rupture has changed. In addition, the test results show that the new method can still provide a more accurate spatial and temporal location even if the focal depth changes obviously.

Testing of the Nepal-Gorkha Earthquake
We integrate the back-projection results of the five arrays to reconstruct the rupture process of the Nepal-Gorkha earthquake using the MCBP method, as shown in Figure 7. In the first half of the rupture, that is 0 ∼ 40 s, the results obtained using our new method are in good agreement with those obtained using the slowness correction back-projection method (Meng et al., 2016) and other back-projection methods (Fan and Shearer, 2015;Liu et al., 2017) because of the main shock correction effect. Compared with the model from GPS and interferometric synthetic aperture radar data (Galetzka et al., 2015), it was found that the location of the high-frequency radiation source released by us was distributed along the fault slip boundary, which just showed that the high-frequency energy was mainly released by the rupture front. The direction of rupture begins to change to the south after the rupture lasts for about 40 s. Considering the strike of the Main Himalayan Thrust (MHT) fault (Lavé and Avouac, 2000), we can infer that the rupture began to propagate to the shallow part of the fault. In the area where the direction of rupture has suddenly changed, Whipple et al. (2016) found steep fault segments. The geological findings of Hubbard et al. (2016) suggest that the northeastern slips reached and ruptured a ramp. Zhang et al. (2017) investigated the 2D dip variations of the Nepal earthquake fault using a new geodetic inversion method and found a dip anomaly in this area. The highfrequency energy radiation location estimated using our MCBP method just bypasses the dip anomalous region, which can be interpreted as the large dip angle anomaly preventing the deep fault from propagating to the southeast. At the later stage of the rupture process, the high-frequency energy is mainly released in the shallow area, but the energy is obviously weakened. It is possible that the shallow part of the fault has a slight slip.
Although the Nepal earthquake was suspended by the steep fault segments, a large number of aftershocks and creep occurred subsequently, which further released the elastic potential energy accumulated in the locking area due to the blocking of the main shock (Elliott et al., 2016). Our results provide new evidence for understanding the fine structure of the Nepal earthquake faults.

DISCUSSION AND CONCLUSION
How to make full use of the information of multiple arrays and integrate their positioning ability has always been a difficult problem in the back-projection of multiple arrays. The direct superposition of the back-projection energy of multiple arrays cannot fundamentally solve this problem. In this study, the error distribution characteristic of back-projection imaging of the array is discussed, and a multi-array back-projection method is developed. The new method makes full use of the azimuth coverage information of each array located in different directions of the source. It improves the positioning accuracy. It solves the problem that the back-projection results of multiple arrays cannot be effectively fused and calibrates the directional systematic error of the single array back-projection.
The seismic ray path changes when the epicenter distance changes. For example, it may pass through different strata or inhomogeneous bodies. Such changes are random and may be discontinuous and non-monotonic, so it is difficult to estimate the change rule with a few large aftershocks.  In this study, we start with the extension array response function and analyze the distribution law of the positioning error of the classical array back-projection. Most of the array response functions are elliptically distributed, which indicates that the errors of the two directions are obviously different. We separate the information with the smaller error from the backprojection results of a single array and design a new multiplearray method.
The epicentral distance error of the array back-projection mainly comes from the error of the 3D velocity structure. The dynamic change of focal depth may be another factor for the error of back-projection results because it is assumed that the focal depth is constant in the back-projection. Xu et al. (2009) calculated the ARF of different depths based on the onedimensional earth structure and found that the influence of depth was very small, but for the back-projection of actual data, we could not completely exclude it. For the onedimensional earth structure model, if the rupture propagates to deep levels, the travel time of the seismic wave from the rupture radiation to each station in the global network decreases synchronously, so the position of the back-projection remains unchanged. However, the increase in depth reflects a lag in the time of rupture. If the three-dimensional earth structure is considered, the change of focal depth may cause seismic rays to pass through different velocity layers in different directions, and the effect of depth may not be ignored. The change of the actual focal depth will eventually be reflected in the location error of the epicenter distance. In addition, when the seismic wave rays recorded by the seismic array are concentrated and the velocity anomaly is encountered, all rays can be considered to have passed through it. The synchronous increase or decrease in travel time is equivalent to the change in the epicenter distance. Therefore, the accuracy of the epicenter distance is mainly affected in the end. The MCBP method eliminates the epicentral distance information obtained from classical back-projection. It would eliminate the influence of the 1D velocity structure model and thus the error caused by the change of focal depth to the greatest extent.
When only two arrays are available, the optimal resolution can be obtained if the two arrays are orthogonal in orientation. This is of guiding significance for seismic monitoring in the target area. For large earthquakes such as the 2004 Sumatra-Andaman earthquake, the subevent jump is not very obvious. However, if we examine the medium-intensity earthquake occurring within the plate, the impact of this sudden jump will be so large that the classical back-projection results can hardly be used as a guide for post-earthquake relief. In particular, higher precision is needed to distinguish the source rupture process of earthquakes of magnitude 6 or so. The proposed method weakens the influence of the velocity model error of the medium and integrates the positioning results of different arrays in physical essence, thus obtaining higher positioning accuracy. This high-precision multi-array back-projection imaging method can be used to study the source rupture process of complex earthquakes. Compared with the back-projection of a single array, our method further improves the accuracy of the results, which is expected to be used to analyze the different frequencies of large earthquakes.
A large number of back-projection imaging studies on the rupture process of the Nepal earthquake show that the existing multi-array back-projection method is effective. The new method we propose is an upgrade of the existing method. The new method strongly relies on the accurate imaging of a single array, so we still recommend developing a deconvolution back-projection method or constructing a perfect array to remove the ARF effect, so as to improve the accuracy of backprojection imaging of a single array as much as possible. If an array with different orientations can give the same result, our method is exactly the same as the result of direct superposition. If the results of each array are inconsistent, our method can extract more accurate information from the epicentral distance and back azimuth for cross-positioning. Therefore, our new method is an upgraded version of the existing multi-array imaging stacking method, which is compatible with the existing methods.
The problem of stacking weight involved in the method of directly stacking multiple arrays is still unavoidable in our method. Due to the difference of aperture, station spacing, and ambient noise, the imaging accuracy of different arrays must be different. In this study, the final results were calculated in the way of equal weight. In the following studies, we will further discuss the issue of weight.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Materials; further inquiries can be directed to the corresponding author.