Spatio-temporal analysis of East Asian seismic zones based on multifractal theory

1 Introduction

Nowadays, humankind faces many environmental problems [1,2,3,4]. Earthquake prediction has always been an internationally recognized worldwide problem [5,6,7,8]. In recent years, with the continuous development of computer science and control theory [9,10], especially nonlinear dynamics theory [11], geological statistics analysis [12] has been rapidly developed. Based on various models and algorithms, researchers have conducted diversified analyses of historical seismic data and obtained many beneficial conclusions [13].

Because seismic source events can be regarded as point events in time or space, seismic events are usually treated as point sets in most seismic models based on nonlinear theory. Since Mantelbrot [14] first introduced the concept of fractal in 1976, fractal geometry, as one of the powerful tools for studying complex phenomena in nature, has been attracting significant attention of scholars worldwide and has been much widely used. In recent years, fractal and multifractal theories have been applied in seismology [15]. These have been fully applied in recognition of precursors of major earthquakes, analysis of aftershock sequence and other fields in relevant seismic regions [16].

In studying the temporal and spatial transmission characteristics of earthquakes, the fractal dimension theory mainly focuses on the fractal structure of the temporal and spatial distribution of earthquakes and their evolution process [17]. As early as in the 1980s, Kagan and Kuopoff used fractal theory to study the spatial distribution of seismic sources and concluded that seismic phenomena are fractal, that is, with scale invariance, at least in a specific space and time range and a certain scale area [18]. However, soon after, the research results of scholars such as Glickman from the Soviet Union, Tadashi Hirayashi from Japan, and Longxing Hirada from Japan all agreed that the spatial and temporal distribution structure of seismic activity was a complex fractal rather than a simple fractal, so the method of multifractal was also introduced into the study of the spatial and temporal characteristics of earthquakes [19,20]. Richa Jain et al. studied the aftershocks of the Channoli earthquake in Garhwal Himalaya in 1999 and found that the correlation dimension of the spatial sequence reached 1.64, which was close to the two-dimensional plane distribution, while the time series became continuous distribution (M ≥ 1) and discrete distribution (M ≥ 1.75) at different magnitudes [21].

With increasing breadth and depth of the application of multifractal theory in seismic-related fields, the relevant models and algorithms are constantly improved and developed in practical applications. Michas et al. [22] used a multifractal approach to study the time dynamics of the recent earthquake activity in the Corinth rift. The results indicated the degree of heterogeneous clustering and correlations acting at all time scales that suggest non-Poissonian behavior. Additionally, the multifractal analysis in different time periods showed that the degree of multifractality exhibits strong variations with time, which are associated with the dynamic evolution of the earthquake activity in the rift and the transition between periods of high and low seismicity. Based on multifractal detrended fluctuation analysis (MF-DFA), Yang et al. [23] calculated the multifractal characteristic parameters of acceleration time series for typical near-fault ground motions and analyzed their correlations with two period parameters (i.e., mean period T _m and characteristic period T _c) and box-counting fractal dimensions. Chamoli and Yadav [24] used the MF-DFA to understand the multifractal behavior of earthquake data. For this purpose, a complete and homogeneous earthquake catalog of the period 1965–2013 with a magnitude of completeness M _w of 4.3 was used. The analysis revealed the presence of multifractal behavior and sharp changes near the occurrence of three earthquakes of magnitude (M _w) >6.6, including the October 2005 Muzaffarabad–Kashmir earthquake.

After a major earthquake, a series of aftershocks generally occur, reflecting the instability of the newly formed geological faults and the discontinuity of energy. In recent years, the fractal theory has also made breakthroughs in this field. Liu et al. [25] applied detrended fluctuation analysis (DFA) and multifractal methods to the time-scaling property analysis of the waiting time sequence of aftershocks (WTSAs) of the Wenchuan Ms8.0 earthquake that occurred on the Longmen mountain tectonic zone of Sichuan Province in China. The results showed that the WTSA is characterized by multifractal scaling.

The propagation characteristics of an earthquake are an interdisciplinary subject. Many research studies combine fractal theory with other theories to study earthquakes. Subhadra et al. [26] characterized the nature of the spatial distribution of heterogeneities in the source area of the 2015 Nepal earthquake based on the seismic b-value and fractal analysis of its aftershocks. The earthquake size distribution of aftershocks gives a b-value of 1.11 ± 0.08, possibly representing the highly heterogeneous and low-stress state of the region. Pasten [27] applied the fractal theory in order to describe the phase transition of earthquakes and had studied the multifractality of the time series of connections on the complex network, and found monofractal behavior for an earthquake with precursors and a multifractal behavior in a region without large earthquake with Mw 6.0. Denisse [28] introduced a new method of characterizing the seismic complex systems using a procedure of transformation from complex networks into time series. It was shown that the multifractality of constructed connectivity time series changes due to the particular geophysics characteristics of the seismic zones – it decreases with the occurrence of large earthquakes – and shows the spatiotemporal organization of these seismic systems.

In this article, the east Asia seismic belt was divided into regions. Then the time series set and the space distribution series set of the five major seismic belts were analyzed by the multi-fractal method. Based on these two sets, we explained the multifractal characteristics of the five regions in seismic time series and seismic spatial distribution series from different aspects. Besides, we elaborated that the multifractal categories of each region are according to the important parameters of their respective relations.

2 Study area and data

Located at the confluence of two major seismic belts, Asia encompasses the eastern part of the Eurasian earthquake zone and the northwestern part of the Pacific Ring of Fire. Due to its complex geological features and high earthquake occurrence frequency, it is highly representative in regions with high earthquake occurrence in the world. Therefore, historical seismic data and modern seismic data of the entire Asian region are of high research value [29]. Based on this, the study area was selected in the range of longitude 80–140°E and latitude 20–60°N, and the seismic data were collected from the SEISMIC catalog of USGS. In terms of time span, the data range was from January 1, 1970, to May 30, 2014, and the minimum magnitude was M ≥ 2.8.

In addition, among the objects studied in this article, the magnitude range was M ≥ 2.3. Rather than discussing only the destructive strong earthquakes (M ≥ 7.0), the Gutenberg–Richter (G–R) relationship can be used to obtain a reasonably complete magnitude catalogue, which makes the research results more objective.

In this article, the magnitude–frequency distribution principle of the G–R relationship was adopted to conduct seismic data integrity analysis [30]:

where M is the earthquake magnitude, N is the number of earthquakes with a magnitude greater than or equal to M in the study area, and a and b are statistical constants.

The steps of seismic directory integrity analysis are as follows:

First, aftershocks of earthquakes (1.0 ≤ M ≤ 9.0) in the study samples (set M ≤ 5.0) will be deleted, as shown in Table 1. A total of 15,138 small earthquakes (M ≤ 5.0) were obtained.

Then, the complete data after the aftershock were deleted according to the step size of 0.1 magnitude compared with Table 2.

After obtaining these data excluding aftershocks, the logarithmic distribution relationship of magnitude frequency of the study samples was plotted, as shown in Figure 1.

Figure 1

Logarithmic distribution of earthquake magnitude frequency in the study area.

The second step was to take different magnitude-lower limits M i ( i = 2.3 − 5.0 ) according to the G–R relation and use the least-square method for fitting to obtain different fitting correlation coefficients.

Through the least-square method [31], the fitting correlation coefficient is given as follows:

The third step was to find the M corresponding to the maximum fitting correlation parameter. When it reaches its peak, M i = M c .

In fact, the whole process is an atlas with a step size of 0.1 from M = 2.3 to M = 5.0. Only representative cases of M = 2.5, M = 3.5 and finally obtained M = 4.5 are listed here, as shown in Figure 2(a)–(c) to illustrate the situation.

Figure 2

Logarithm of magnitude frequency in East Asia and its fitting G–R figure. (a) M = 2.5, (b) M = 3.5, and (c) M = 4.5.

As can be seen from Figure 2(a), when M = 2.5, the fitting correlation coefficient R = 0.763, which is low and the fitting linear relationship is poor. When M = 3.5, as shown in Figure 2(b), the fitting correlation coefficient has reached a relatively high value. At this point, R = 0.836, and the fitting linear relationship is much better than before. When M = 4.5, as shown in Figure 2(c), the fitting correlation coefficient R = 0.997 has reached a very high level, indicating that the fitting linear relationship is very good.

In the further experiment, when M ≥ 4.6, the value began to decline significantly. According to the theory, this is because with an increase in the value of M c , the data of smaller magnitude–frequency relation distribution with good linearity are gradually discarded, which reduces the reliability of the fitting result. Therefore, it can be determined that when M = 4.5, M i is M c . The change in the process of fitting correlation coefficient R with the value of M in the whole process is shown in Figure 3.

Figure 3

Diagram of fitting correlation coefficient R changing with the value of M.

Therefore, excluding the seismic data of M < 4.5 from all previous data, this is our final research data.

Due to the geographical scope of this study [32], we can’t use a single fractal and fractal dimension of parameters to describe. Even if getting ideal results under these methods, the concluded experimental results do not have a big reference value on regional earthquakes due to the combination of various regions. In addition, as the topography and landforms in East Asia are very complex, the whole sample will be regionally divided according to different geological structures.

Figure 4 shows the southeast corner of China, and each black dot in the image represents seismic data. We first divided China’s seismic belt into two large areas according to the blue line in Figure 4. To the left of the blue line is part one, and to the right of the blue line is part two. The terrain and geomorphology of the region where the first part is located are relatively complex, which cannot be divided into a separate whole but needs to be divided further.

Figure 4

The regional division diagram of the research object.

The division of the first part is mainly based on the seismic zone of mainland China. Seismic research in this article cannot easily divide all the traditional seismic zones for all the data included in the study, but it roughly divides according to the alignment of these areas, respectively: the seismic zone as the center, the center into the surrounding area formed by the division into blocks, each block in the separate seismic data including the research object, the piece of non-interference in each other’s data in different areas. The general division of the research area is shown in several sections in Figure 5. The figure shows the entire area of China. The four ends of China’s territory are as follows: the easternmost is at the intersection of the Central Line of the main waterway of Heilongjiang and Wusuli River (135°2′30″E), the westernmost is near the Pamir Plateau (73°29′59.79″E), the northernmost point is on the center line of heilongjiang Main Waterway north of Mohe (53°33′N, 124°20′E), and the southernmost is at the Ground ansha (3°31′00′N′, 112°17′09″E). Lidi Ansha or Lydi Shoal is the southernmost part of the de facto Chinese territory. At the lower right corner of the map is a map of South China Sea territory. If these are all represented in a large map, it will lead to uncoordinated proportions or a too large picture, so these are represented separately [33].

Figure 5

Schematic diagram of major earthquake zones in mainland China.

The first part is divided into four seismic zones, that is, the south–north earthquake zone of China, the earthquake zone of Xinjiang and its surrounding areas, the earthquake zones of the Qinghai–Tibet Plateau and The Himalayas, and the earthquake zone of north China. At the same time, the second part as a whole is an earthquake zone, known as the northwest Pacific Rim earthquake region.

3 Methods

The fractal theory has an important basic characteristic, namely, self-similarity [34]. Self-similarity refers to symmetry across different scales (including time scale and space scale).

Self-similarity or scale invariance can be mathematically expressed as follows:

A simple function that satisfies this property is a power function [32]:

Thus, it can be obtained that [35]:

For regular fractals, the scale-extended multiple λ is the intermittent value.

In seismology, the empirical relations to describe the main characteristics of seismic activity, such as Omori’s formula that the aftershock attenuates with time and the g–R relation, are distributed in the form of a power function. When Kagan and Kuopoff studied the local spatial distribution [35], they found that the number of earthquakes per unit volume, N, also has the following relationship with the focal distance r over a relatively wide scale:

Therefore, there is scale invariance between many parameters in seismology.

In determining simple fractal dimensions of irregular fractal, the box counting method is often used. The multifractal theory considers the slightly different physical quantities within the box. It is normalized by a statistical method and described by multifractal spectrum [36]. In this way, it is possible to restore the distortion of the final result caused by ignoring the differences within each box.

Renyi generalized that fractal dimension is the most common way to describe multifractals. It defines a q-degree probability moment from the perspective of information theory [15]:

When ε → 0 , there is

where D q is the generalized dimension. Its definition is divided into two parts according to different q values.

When q = 1

When q ≠ 1

The selection of power parameter q (hereafter referred to as power parameter) is very important. The parameter defined in theory can be taken over the whole real number. The larger the range, the better, it can fully show the fractal dimension change. However, in practical application, with the increase of numerical value, the amount of calculation will also increase in the form of power, resulting in the overflow of calculation results and the failure of drawing. If the value range of the parameter is too small, the complete multifractal spectrum curve cannot be obtained, resulting in the loss of a significant amount of necessary information, and the multi-fractal characteristics of the object cannot be completely studied. Therefore, the selection of q needs to be determined by combining experience values and experiments.

Another set of methods to describe multifractals is the multifractal spectrum. The research object was divided into N regions S i according to scale ε i , and the corresponding probability was P i ( i = 1 , 2 , … , N ) . When S i is small enough, P i = ε α i can be considered. Here α i is called the scaling exponent, or singular strength. In many α i , the number of regions N i ( ε ) with the same value α is also related to scale ε i , that is [37]

Here, f(α) is the component of alpha in the total distribution, and it is the continuous derivative of α . The multifractal spectrum, f(α) spectrum, is constructed.

Through Legendre transformation, the following can be obtained:

Let τ ( q ) = q a − f ( α ) , then

It follows that D q − q and f ( α ) − α are equivalent and can be solved interchangeably.

4 Results

In the study of multifractal, the sample size cannot be ignored. Many studies have shown that when the sample size is less than 200 times, the multifractal dimension presents a pseudo-multifractal characteristic with the variation in parameters, while when the sample size is much greater than 200 times, the calculated results can objectively reflect the properties of the samples [38]. By comparing the sample data of various regions in Table 3, it can be predicted that the amount of seismic data in the north China earthquake region is less than 200, so the seismic data set in this region cannot objectively reflect the multifractal structure either from the time series or from the spatial distribution series and may even show the pseudo-multifractal characteristics.

During the study of seismic time series in this article, according to the selected seismic data and combining the empirical value and the results of multiple attempts, the final determined |q| is at most 10 and increments as ∆q = 0.1 in the range [−10,10]. The multifractal spectral curve of f(α) can approximately reflect the changing characteristics within the q ∈ (−∞, +∞) range.

In the specific model calculation, we first calculated the double logarithm ratio under different q values of the seismic time series obtained by the algorithm of counting box dimension. Figure 6 shows the logarithmic atlas of the time series of each earthquake region. Because when the power parameter is greater than 4, the linear relationship of the log–log graph in some areas becomes worse, for the convenience of comparison, we temporarily take q ∈ [−4,4] in the graph, and draw the ln X(q)–ln L graph with a step length of 0.5, where X(q) is the partition function, which is the sum of the q-power of the probability distribution set after the normalization of points in the time series. Each power parameter corresponds to a logarithmic curve.

Figure 6

ln X ( q ) − ln L diagram of time series of each seismic region. (a) North-South, (b) Xinjiang, (c) Qinghai-Tibet, (d) North China, and (e) Pacific Rim.

It can be seen from Figure 6 that no matter the power parameter is less than or greater than 0, all the ln X ( q ) curves in L intervals show a good linear relationship with ln L , indicating that the time series of each region belongs to a multifractal. With a section of time series ln X ( q ) − ln L diagram, compared to corresponding spatial distribution sequence diagram, when the power parameter is greater than 4 linear gradually becomes poor. This is due to the spatial distribution of sequence data grid partition and L is the minimum limit. When the power parameter is greater than 4, the amount of sample data reduces sharply, so the linear relationship between decreases significantly that the multifractal characteristics cannot fully exhibit, so the double logarithm figure appeared certain fluctuation.

In the figure, when the power parameter value is > 0, ln X ( q ) can change from −10 with ln L and then basically remains as a straight line. However, when the power and value are <0, ln X ( q ) can only start from ln L = 4 values; this is because in the process of calculation, the value of ln L is too small and the timeline e is segmented into so many segments that each segment of the seismic data of probability will change especially small, small to the computer here instead of relying on their tail. Therefore, ln X ( q ) take the outliers. In fact, if the calculation accuracy can ensure that this mantissa is not discarded here, when the power parameter is <0, when the value of ln L is small, ln X ( q ) will be much larger than 25. With the increase of ln L , ln X ( q ) will rapidly decrease. This is a process of abnormal decline. On the other hand, not all L ranges meet the scale invariance.

After obtaining the logarithmic map of different power parameters in each region and verifying the multifractal characteristics of the earthquake time series, the slope of the ln X ( q ) − ln L curve can be used to calculate τ ( q ) . The relationship between τ ( q ) and the partition function is as follows [39]:

If the object of study is an irregular multifractal, the slope of τ ( q ) is piecewise; that is, the change of τ ( q ) with q is nonlinear. The stronger the nonlinearity of τ ( q ) is, the stronger the multifractal strength is. Due to the limited sample space, the calculation method here is to draw the logarithmic graph of ln X ( q ) − ln L in common use according to the power parameter value and then use the invariable linear regression method to fit and calculate a series of data, namely,

Finally, the slope of the line can be obtained as follows:

Then, according to the given power parameter value, we perform unitary linear regression fitting for each curve, calculate the value of τ ( q ) , and draw its corresponding power parameter distribution diagram. As shown in Figure 7, the slopes of all the curves are obviously divided into two sections, indicating that the generalized energy time series of the earthquake does have multifractal characteristics, and the difference of the slopes represents the strength of the degree of multifractal.

Figure 7

The comparison of time series τ ( q ) − q in each earthquake region.

After the ln X ( q ) − ln L relation and the τ ( q ) − q relation of time series in each earthquake region are obtained, two sets of multifractal spectra can be obtained to describe the multifractal characteristics of each earthquake time series according to the multifractal theory [15].

First, let us focus on the D q − q relation. According to the multifractal theory, from formulas (9)–(11), we found that when the power parameter is > 0, the larger P i plays a major role, and D q mainly describes the structural complexity of the cluster area of the earthquake time series. When the power parameter is <0, the smaller P i plays a major role, and D q mainly describes the structural complexity of the sparse region of the seismic time series.

The D q − q relation of the time series of the five seismic regions is shown in Figure 8. With the gradual dispersion of the samples in the sparse region, the D q − q curve of the cluster in the dense region increases, D − ∞ increases, and D + ∞ decreases.

Figure 8

Time series D q − q diagram of each earthquake region.

It can be seen that D + ∞ is stable over D 0 and D q of the negative power parameter. Many studies have shown that d infinity decreases significantly before a major earthquake with the change of time and tectonic energy. However, the research object of this article is relatively broad in geographical scope and the singularity of small local scope is averaged out, so this phenomenon is not very obvious. According to the dynamic characteristics of the D q − q curve, the time series distribution of the general seismic background activity D q − q shows slow deformation. However, if the earthquake in this area is more complex, the D q − q relation curve will show a steep deformation, and D q decreases sharply with an increase of q. Table 4 lists several important parameters of the D q − q relationship curve.

Among them, ΔD can reflect the size of the multi-fractal complexity of a region. It can be seen that the values for north China, Qinghai-Tibet, and Himalayan earthquake regions are relatively high, while the northwest seismic region around the Pacific Ocean has the lowest value.

As mentioned before, due to the limitation of sample size, it is impossible for the north China seismic region to reflect the accurate multifractal characteristics in the algorithm of this article. Therefore, the high value of ΔD in north China is actually a pseudo-multifractal property. The rest of the region as the data samples are suitable and ensure the accuracy of the results. The high value of ΔD in the Tibetan and the Himalayan seismic area may be due to the uneven energy release cased by the regional tectonic belt under various stress, resulting in the higher multifractal complexity. The low values in the northwestern Pacific Rim earthquake zone may be due to a more uniform release of seismic energy in the region.

Non-uniform complex fractals need to be characterized by this set of dimensions { D q }. These can also be characterized by another language scale index spectrum, that is, the f ( α ) − α relationship. We can associate D q with f ( α ) − α by using a Legendre transform. Figure 9 shows the f ( α ) − α relationship for several seismic areas.

Figure 9

Time series f ( α ) − α diagram for the northwestern Pacific Rim earthquake region.

It can be seen that although the D q − q diagram can reflect the difference between the dense and sparse multifractal regions to some extent, the effect is not as obvious as the f ( α ) − α diagram. According to the morphological classification of the f ( α ) − α graph, the curve can be divided into three types. The first is sparse, with vertices skewed to the left. The second is intensive, with vertices biased to the right. The third is compound, where the vertex is centered, and the curve is bell-shaped.

By comparing Figure 9, it can be seen that the seismic time series of the earthquake regions in Xinjiang and surrounding areas, the Qinghai–Tibet Plateau, and the Himalayan earthquake regions have the obvious characteristics of sparse multifractals. The seismic time series in the north China seismic region is slightly inclined to the intensive multifractal characteristics. The seismic time series of the north and south seismic regions of China and the northwest Pacific Rim are complex.

Unlike time series, seismic spatial distribution data need to be prepared for reserialization. The region of the grid is determined. In the processed data of the partition, the longitude and latitude information of each earthquake is included into a set, respectively, and the maximum and minimum values of the longitude and latitude of the set are found, thus determining the boundary of the rectangular region. Then, the rectangle is meshed in the rectangular area with a side length of 0.2°, so that each point of the seismic distribution falls into the grid. Then, the points of the grid are counted, and the set is serialized according to the sorting of the grid. After obtaining the sequence, take q ∈ [−4,4] and draw the ln X(q)–ln L diagram with a step length of 0.5.

It can be seen from Figure 10 that no matter the power parameter is less than or greater than 0, all the curves ln X ( q ) in L intervals show a good linear relationship with ln L , indicating that the space series of each region belongs to a multifractal. With a section of time series ln ln X ( q ) − ln L diagram, compared to corresponding spatial distribution sequence diagram, when the power parameter is greater than 4 linear gradually becomes poor. This is due to the spatial distribution of sequence data grid partition and L is the minimum limit. When the power parameter is greater than 4, the amount of sample data reduces sharply, so the linear relationship between decreases significantly that the multifractal characteristics cannot fully exhibit, so the double logarithm figure appeared certain fluctuation.

Figure 10

The spatial distribution sequence ln X ( q ) − ln L diagram of each seismic region. (a) North-South, (b) Xinjiang, (c) Qinghai-Tibet, (d) North China, and (e) Pacific Rim Tibetan.

In Figure 10, when the power parameter value is >0, ln X ( q ) can change from −10 with ln L and basically remains as a straight line. However, when the power and value are <0, the same happens as the time series analysis since not all L ranges meet the scale invariance, and ln X ( q ) can only start from ln L = 4 values. After obtaining the logarithmic map of different power parameters in each region and verifying the multi-fractal characteristics of the seismic spatial distribution sequence, the slope of the ln X ( q ) − ln L curve can be used also to calculate τ ( q ) .

As shown in Figure 11, the slopes of all the curves are obviously divided into two sections, indicating that the generalized energy time series of the earthquake does have multifractal characteristics, and the difference in the slopes represents the strength of the degree of multifractality. Figure 11 shows the comparison of the spatial distribution sequence τ ( q ) − q in each earthquake region. After obtaining the relationship between ln X ( q ) − ln L and τ ( q ) − q of the spatial distribution sequences of each earthquake region, two sets of multifractal spectra can be obtained to describe the multifractal characteristics of the spatial distribution sequences of each earthquake according to the multi-fractal theory.

Figure 11

The comparison of spatial distribution sequence τ ( q ) − q in each earthquake region.

Thus, the spatial distribution sequence D q − q of the five seismic regions is obtained, as shown in Figure 12. With the gradual dispersion of the samples in the sparse region, the D q − q curve of the cluster in the dense region increases, D − ∞ increases and D + ∞ decreases.

Figure 12

Spatial distribution sequence D q − q diagram of each seismic region.

According to the dynamic characteristics of the D q − q curve, the spatial distribution sequence of general seismic background activity D q − q shows slow deformation, and D q decreases slowly with an increase of q . However, if the earthquake in this area is more complex, the D q − q relation curve will show a steep deformation, and D q decreases sharply with an increase of q . Table 5 lists several important parameters of the D q − q relationship curve, where the magnitude of ΔD can reflect the magnitude of the multifractal complexity of a region. Similarly, due to the limitation of the sample size, the multifractal spectrum of the north China earthquake region has no reference significance. However, there is no limit of minimum time interval in the analysis of the Northwestern Pacific Rim earthquake zone, which can better reflect the multi-fractal characteristics of spatial distribution.

It can be seen that the values of the Tibetan and Himalayan seismic regions are relatively high, which is consistent with the time series, indicating that the Tibetan and Himalayan tectonic belts are not uniformly released due to various stresses, resulting in a high multifractal complexity. The values of the Northwestern Pacific Rim earthquake zone are also high, which indicates that the energy release of this tectonic belt is more complex, which is in accordance with its geological tectonic background. The energy release in Xinjiang and its surrounding earthquake areas and the earthquake areas in the north and south of China is much more uniform than the previous two.

Figure 13 shows the f ( α ) − α diagrams for several seismic areas. By comparing Figure 13, the characteristics of regional intensive multifractals in north and south China earthquakes can be seen. The peak of seismic spatial distribution in the Qinghai–Tibet Plateau and the Himalayan region is left-biased, which is a typical sparse multifractal characteristic. The peaks of the seismic region of Xinjiang and its surrounding area with the Northwestern Pacific Rim seismic region are bell-shaped, which belongs to complex characteristics.

Figure 13

f ( α ) − α diagram of the spatial distribution sequence of the northwest Pacific Rim earthquakes.

5 Discussion

In this article, the seismic catalog of seismic data in the longitude range of 80–140°E and latitude range of 20–60°N from 1970 to 2014 in the Chinese mainland is used as the basic research data. After data preprocessing, the study area is divided into separate research areas. The multifractal model of time series and spatial distribution series of each seismic region is established, and the logarithmic relationship graph of each seismic region with a power parameter of [−4,4] and a step size of 0.5 is shown. After that, the τ ( q ) − q relation, D q − q relation, and f ( α ) − α relation are further established to illustrate the multifractal characteristics of five regions in earthquake time series and seismic spatial distribution series from different aspects, and the multi-fractal categories of each region are elaborated in detail according to the important parameters of each relation.

Because the multifractal analysis based on the number box method has a minimum limit on the size of the sample size, a collection of <200 cannot accurately reflect the regional multifractal properties, so the sample size is only 169 out of the north China earthquake area, the results show that the area both in time series and spatial sequence shows the pseudo multifractal features, also certain degree to verify the theory need the support of large sample size. The following content focuses on the multi-fractal analysis of the remaining regions from the two perspectives of the relationship between D q − q and f ( α ) − α . The results show that from the perspective of the time series D q − q spectrum, the high values in the Tibetan and Himalayan seismic regions may be due to the uneven energy release caused by various stresses in the tectonic belt of these regions, resulting in the high multifractal complexity. However, the multifractal complexity of Xinjiang and its surrounding earthquake regions as well as the north–south earthquake regions of China is obviously weaker than that of the Qinghai–Tibet and Himalayan earthquake regions, which indicates that the activities of these two tectonic belts are not as intense as those of the Qinghai–Tibet and Himalayan earthquake regions. From the perspective of the relationship between time series f ( α ) − α , the seismic time series between Xinjiang and its surrounding regions and the Tibetan Plateau and Himalayan regions are characterized by sparse multifractals. The seismic time series in the north–south seismic region of China shows the complex characteristics. From the perspective of the D q − q spectrum of the spatial distribution sequence, the values of the Tibetan and Himalayan seismic regions are higher and consistent with the time series, which indicates that the Tibetan and Himalayan tectonic belts are affected by various stresses, resulting in uneven energy release, resulting in a high multifractal complexity. The values of the Northwestern Pacific Rim seismic region are also high, which indicates that the energy release of this tectonic belt is more complex, and this result is consistent with the geological and tectonic background. The energy release in Xinjiang and its surrounding earthquake area is much more uniform than that in the north–south earthquake area of China. From the point of view of the relationship of spatial distribution sequence f ( α ) − α , the regional intensive multifractal characteristics of earthquakes in North and South China are studied. The seismic spatial distribution vertices in the Tibetan Plateau and Himalayan earthquake regions are left-skewed, which belongs to the typical sparse multi-fractal characteristics. The seismic spatial distribution vertices in Xinjiang and its surrounding earthquake area and the Northwestern Pacific Rim seismic region are bell-shaped, belonging to the complex characteristics.

In the process of analysis, some defects of this study have also been fully exposed. First, the study of the earthquake region in north China cannot accurately reflect the characteristics of seismic space–time propagation due to the limitation of the sample size for multifractal analysis, so this region must be discussed by other methods. There are three methods to calculate a generalized fractal D spectrum with multifractal data, that is, number box method, fixed radius method, and fixed mass method. For the region, this study has abandoned the number box method. Besides, the fixed radius method for the sparse earthquake area has a big statistical error, which cannot be used. The methods of fixed quality and fixed radius asymptotic are equivalent but present the opposite statistical behavior; therefore, the area is suitable for further use for fixed quality standard of the multifractal spectrum analysis. The establishment of the model is very important, and the final calculation result depends on the specific calculation process in the model to a large extent. In this article, the selection of the minimum interval of the time series, the division of the maximum rectangular range of the spatial distribution sequence, the selection of the minimum grid and the final serialized grid order have a great impact on the calculation results, so the calculation accuracy cannot be determined finally. However, from the domestic and foreign research in this field, it is also a huge obstacle to the development of this field, so whether these fractal dimensions can more accurately describe the spatial and temporal propagation characteristics of earthquakes in a region has yet to be determined by the fractal dimension theory itself.

6 Conclusion

First, this article analyzes the catalog integrity of the original samples according to the G–R relationship. The fractal analysis has a high requirement on the size of the sample. In a fractal model, the required amount of data shall not be less than 200. In the research of earthquakes, especially the more remote area, due to the lack of historical data, many researchers do not have complete historical data, resulting in the final analysis result lost its reference significance. As the research object in this article has a large time span and a wide geographical scope of the research area, the data are relatively rich. In this case, the integrity of the original data is analyzed in this article and the minimum magnitude of the research object is determined.

Second, because the geographical scope of the object studied in this article is very wide and its geological background is very complex, if the unified fractal dimension is used as a parameter to describe the space–time characteristics of the whole object, the result will be distorted due to averaging out the local singularity. In order to more accurately describe the fractal dimension characteristics of each region, this article divides the research region into five research regions according to different geological structures and seismic zones.

On this basis, this article analyzes and calculates each area according to the theory of multifractal theory and discusses the final results according to the changes of important parameters. The multifractal model of time series and spatial distribution series of each seismic region is established respectively, and the logarithmic relation diagram of each seismic region with power parameter of [−4,4] and a step size of 0.5 is displayed. The τ ( q ) − q relation, D q − q relation, and f ( α ) − α relation are further established. The characteristics multifractal in seismic time series and seismic spatial distribution series of five regions are explained from different aspects, and the multifractal categories of each region are described in detail according to the important parameters of their respective relations.