JP6891001B2 - Seismic intensity distribution estimation method and estimation system - Google Patents

Description

The present invention relates to a method and an estimation system for estimating the seismic intensity distribution, which makes it possible to estimate the distribution of the seismic intensity over a wide area at various timings with high accuracy after the occurrence of an earthquake.

Accurately predicting the area distribution of seismic intensity over a wide area when an earthquake occurs is to understand structures and human damage, establish rescue and evacuation postures, and formulate recovery plans. Is important for.

Therefore, in recent years, based on the observation data of earthquakes that have occurred in the past, the wide-area distribution of seismic intensity such as maximum acceleration and velocity is estimated in real time after the occurrence of the earthquake or immediately after the occurrence of the earthquake. Various systems are being developed.

In the conventional estimation system of this kind of seismic intensity distribution, as generally disclosed in Non-Patent Document 1 below, a distance attenuation formula obtained from regression analysis based on observation data of seismic motions that have occurred in the past. By estimating the seismic intensity on hard ground and multiplying it by a coefficient that evaluates the easiness of shaking of the surface layer, a method of estimating the seismic intensity analysis on the ground surface is adopted. Most of the damage estimations carried out by local governments use this method to estimate the seismic intensity.

However, the distance attenuation equation used in the estimation system is basically a function of magnitude and distance, as shown in the equation below.
Distance attenuation formula: logY = aM-logX-bX + c (1)
Here, Y is the seismic intensity such as the speed of shaking at the estimated location, M is the magnitude at the epicenter position, and X is the distance to the epicenter.

Therefore, in the prediction formula in which the coefficients a, b, and c are determined by regression analysis using the above formula (1), the distribution of seismic intensity is concentric around the epicenter position, whereas in the past, for example. In a study on the seismic intensity distribution of earthquakes that occurred in the Tohoku region, it was pointed out that the distribution shape does not become concentric around the epicenter position, such as the seismic intensity distribution extending from north to south (abnormal seismic area). There is.

That is, Fig. 4 shows the distribution of the maximum velocity obtained by the seismic observation network of the National Research Institute for Earth Science and Disaster Prevention during the M5.6 earthquake that occurred in Ibaraki Prefecture at 18:56 on March 19, 2011. Is. It can be seen that even in this earthquake, the distribution of the maximum velocity values is not concentric around the epicenter position, but extends north and south.

On the other hand, FIGS. 5 to 7 show the results of estimating the distribution of seismic intensity for the earthquake shown in FIG. 4 by the above-mentioned conventional method. In this method, the maximum velocity distribution on the engineering basis is first obtained using the distance attenuation formula disclosed in Non-Patent Document 1 (Fig. 5), and then the results are published in the National Research Institute for Earth Science and Disaster Prevention. The maximum velocity distribution on the ground surface was estimated by multiplying the ground amplification factor (Fig. 6) from the target base to the ground surface (Fig. 7).

By comparing FIG. 7 showing the estimated results with FIG. 4 described above, it can be seen that the distribution of the seismic intensity extending from north to south cannot be appropriately evaluated by the above-mentioned conventional method.

Therefore, the present inventors have conducted intensive studies to improve the prediction accuracy of the ground motion intensity distribution using the above distance attenuation formula, and have focused on the fact that there is a mutual relationship between the transmission of ground motion between observation points, and the distance. When determining the coefficient in the decay equation, an evaluation is made in consideration of interregional factors in fields such as economics and statistics, which are proposed in the following Patent Documents 2 and the like and are found in the following Non-Patent Documents 3 and 4, for example. In order to obtain the knowledge that the above-mentioned estimation of the seismic intensity distribution in the abnormal seismic region can be estimated with high accuracy by weighting between the observation points by utilizing the geographical weighted analysis analysis used in the case. I arrived.

Japanese Unexamined Patent Publication No. 2011-96004 Tsukasa, Midorikawa: Maximum acceleration / maximum velocity distance attenuation formula considering fault type and ground conditions, Architectural Institute of Japan Structural Papers, 523, pp.63-70, 1999 A. Stewart Fotheringham, et al .: Geographically Weighted Regression the analysis of spatially varying relationships, JOHN WILLEY & SONS, LTD, 2002 Research on Neighborhood Externalities of Kyomachiya Accumulation, JSCE Proceedings D, Vol.62 No.2 p227-238, 2006.6 Understanding the seasonal and spatial characteristics of the thermal environment mitigation effect of green areas using the geographically weighted regression method, Japan City Planning Association, City Planning Report No.9 p93-97, 2010.11

The present invention has been made based on the above findings, and an object of the present invention is to provide a method and an estimation system for seismic intensity distribution that can estimate the distribution of seismic intensity over a wide area with high accuracy. It is a thing.

In order to solve the above problem, the method for estimating the seismic intensity distribution according to claim 1 is based on the magnitude of past earthquakes, the intensity of seismic motion observed at a plurality of observation points, and the spatial coordinates of the observation points. A prediction formula for seismic intensity is constructed by obtaining the regression coefficient in a regression model that uses geographically weighted regression analysis for the distance attenuation formula, and the distribution of the seismic strength is estimated using the prediction formula when a new earthquake occurs. It is characterized by doing.

Further, the seismic intensity distribution estimation system according to claim 2 is geographically weighted by a distance attenuation formula based on the magnitude of past earthquakes, the intensity of seismic motion observed at a plurality of observation points, and the spatial coordinates of the observation points. A first calculation means for obtaining a prediction formula for seismic intensity by obtaining a regression coefficient in a regression model using regression analysis, a storage means for storing the above prediction formula obtained by this first calculation means, and a new one. Distribution of the seismic intensity of the earthquake by substituting the receiving means for receiving the magnitude and the epicenter position when an earthquake occurs and the distance from the magnitude and the epicenter position received by this receiving means into the above prediction formula. It is characterized by including a second calculation means for estimating the above and a display means for displaying the distribution of the seismic intensity estimated by the second calculation means.

Further, the invention according to claim 3 is the invention according to claim 2, wherein the storage means includes the magnitude of the past earthquake, the strength of the seismic motion observed at a plurality of observation points, and the observation point. The spatial coordinates, the magnitude of the new earthquake received by the receiving means, the strength of the seismic motion observed at the observation point, and the spatial coordinates of the observation point are stored, and the first calculation means is The regression coefficient is obtained from the magnitude of the past earthquake and the new earthquake, the strength of the seismic motion observed at the plurality of observation points, and the spatial coordinates of the observation point, and is stored in the storage means. It is characterized by updating the above prediction formula.

According to the invention according to any one of claims 1 to 3, when obtaining the coefficient of the distance attenuation equation, weighting between observation points is performed using geographic weight analysis analysis, and the optimum regression coefficient at each observation point. By setting, it is possible to estimate the seismic intensity distribution in the above-mentioned abnormal seismic region with high accuracy.

At this time, in particular, according to the invention of claim 3, the magnitude of past earthquakes, the intensity of seismic motion observed at a plurality of observation points, and the spatial coordinates of the observation points are sequentially adjusted to the magnitude at the time of a new earthquake. By updating the above regression coefficient by adding the strength of the seismic motion observed at the observation point and the spatial coordinates of the observation point, the accuracy of the prediction formula can be improved over time.

It is a schematic block diagram which shows one Embodiment of the estimation system of this invention. It is a graph which shows the change of the value of a distance and a space weighting matrix when the bandwidth which determines a weighting range is 50km when space weighting is performed using a bi-square type function in this embodiment. It is a distribution map of the maximum velocity which shows the result of the seismic intensity distribution estimated by the embodiment of this invention. It is a figure which shows the distribution of the maximum velocity obtained by the seismic observation network of the National Research Institute for Earth Science and Disaster Prevention at the time of the M5.6 earthquake that occurred in Ibaraki prefecture at 18:56 on March 19, 2011. It is a distribution map of the result of finding the maximum velocity on the engineering basis for the earthquake of FIG. 4 using the conventional distance attenuation formula. It is a figure which shows the ground amplification factor from the engineering foundation to the earth surface published in the National Research Institute for Earth Science and Disaster Prevention. It is a figure which shows the maximum velocity distribution estimated by multiplying the result of FIG. 5 by the ground amplification factor of FIG.

FIG. 1 shows an embodiment of the seismic intensity distribution estimation system according to the present invention.
In this estimation system, a receiving device (receiving means) 1 for receiving information on the ground motion intensity observed at the observation point i (i = 1, 2, ..., N) at the time of an earthquake occurs, and the receiving means 1 From a connected general-purpose computer (for example, a PC) 2, a storage device (storage means) 3 such as a hard disk connected to the PC 2, and a monitor (display means) 4 that displays the result obtained by calculation by the PC 2. It is roughly configured.

In this PC 2, a RAM, a storage device 3, an input device such as a keyboard or a mouse, and a monitor 4 for displaying input / output data are connected to a CPU (main control unit) that controls the entire system via an input / output control unit. It is well known.

Then, in the storage device 3, the magnitude M of the past earthquake, the strength of the seismic motion observed at the plurality of observation points i (maximum speed in this embodiment), and the spatial coordinates of the observation points i (u _i , v _i ). Specifically, a database of latitude (u _i ) and longitude (v _i ), a prediction formula to be described later calculated using these magnitude M, seismic intensity and spatial coordinates (u _i , v _{i), and the relevant} A first calculation program (first calculation means) for calculating a prediction formula, and a second calculation program (second calculation) for estimating the seismic intensity distribution when a new earthquake occurs from the prediction formula. Means) is stored.

Here, the first arithmetic program is for obtaining the prediction formula of the seismic intensity by obtaining the regression coefficient in the regression model using the geographical weighted regression analysis for the distance attenuation formula. To explain this concretely, first, the distance attenuation formula of the seismic intensity is expressed by a general formula as follows (referred to as a global model).

Here, β ₀ is the regression coefficient that is the constant term (intercept) _{at the observation point i, β k} is the regression coefficient of the kth explanatory variable at the observation point i, x _ik is the explanatory variable (magnitude and source distance), y _i Is the objective variable (maximum velocity value), and ε _i is the regression error.

Then, in the present embodiment, by using the geographical weighted regression analysis, it is allowed that the regression coefficients β ₀ and β _k are different for each observation point i, and the functional type thereof is expressed by a general formula. It becomes like the following formula (2).

Here, (u _i , v _i ) represents the spatial coordinates (latitude and longitude) of the observation point i as described above. In this geographically weighted regression analysis, it is assumed that the seismic intensity y _i of the observation point i is more affected by the nearby observation points than the distant observation points, and the neighborhood observation information of the observation point i is partially sampled. The regression coefficient is estimated by regression analysis using the local weighting.

Then, when the above general formula is applied to the above-mentioned distance attenuation formula (1),
logY _i = aM−logX−β (u _i , v _i ) X ＋ β ₀ (u _i , v _i ) ＋ ε _i (3)
Can be represented by.

Here, X indicates the above explanatory variable, Y indicates a matrix of observed values, and T indicates a transposed matrix. W (u _i , v _i ) is expressed by the following equation as a spatial weight matrix.

The global model shown at the beginning is solved by setting all the elements _{of w i 1} to w _in of this spatial weight matrix to 1. On the other hand, w _ij of the above spatial weight matrix means how much weight is given to the observed value at the j point with respect to the observed point i. Various methods have been proposed for setting this weight value, but the simplest one is 1.0 for the distance between the observation point i and the observation point j that is d or less, and the one that is farther than d. There is a method of setting it to 0.0.

In the present embodiment, the value of the spatial weight matrix is determined by using the bi-square type function shown in the following equation.

In this bi-square function, b (bandwidth) that determines the weighting range is used as a variable. For example, when b = 50 km, the distance dij between the observation point i and the observation point j, The relationship with the value of the spatial weighting matrix w _ij is shown in FIG.

_{Then, the first arithmetic program inputs the magnitude M of the past earthquake and the intensity (maximum velocity value) y i} of the seismic motion observed at the plurality of observation points i into the above equation (3), and also inputs the observation point. Enter the distance X obtained from the spatial coordinates of i and the epicenter position, and select a value _{preset as the bandwidth b in equation (4) above or a value that minimizes the regression error ε i obtained by iterative calculation.} As a result, the prediction formulas in which the coefficients a and the regression coefficients β ₀ and β at each observation point i are determined are obtained and stored in the storage device 3.

In addition, the first arithmetic program is the magnitude M of the earthquake, the intensity of the seismic motion (maximum velocity value) y _i observed at multiple observation points i, and the space of the observation points i after a new earthquake occurs. Coordinates (u _i , v _i ) are added to the above database, and from these data observed at a plurality of observation points i at the time of the past earthquake and the above new earthquake, similarly returned by the above equation (3). It also has a function of obtaining the coefficients β ₀ , β, etc., and updating the prediction formula stored in the storage means 3.

In addition, the second arithmetic program provides the seismic intensity (magnitude) at a plurality of observation points i received by the receiving device 1 and the spatial coordinates (u _i , v _{i) of the observation points i when a new earthquake occurs.} ) Is input to the prediction formula stored in the storage device 3 to estimate the seismic intensity (maximum velocity value) at the observation point i.

Thus, in the method of estimating the ground motion intensity distribution of the present embodiment, the spatial coordinates (u _i strength y _i and the observation point i of ground motion observed in the magnitude M and a plurality of observation points i of past earthquakes , V _i ) and the seismic intensity by obtaining the _{regression coefficients β 0} , β, etc. in the regression model using the geographical weighted regression analysis in the distance attenuation equation using the above equation (3) by the first arithmetic program. Is constructed, and when a new earthquake occurs, the distribution of the seismic intensity is estimated using the prediction formula by the second arithmetic program.

Incidentally, FIG. 3 shows the results of trial calculation of the distribution of seismic intensity using the above equation (3) for the earthquake shown in FIG. However, in this trial calculation, since only one earthquake is used, the magnitude M coefficient a is not obtained. In addition, the bandwidth b of the bi-square function for weighting was b = 44 km by searching for the one with the smallest _{regression error ε i in the program.}

As can be seen in FIGS. 3 and 4, according to the estimation method using the seismic intensity distribution estimation system having the above configuration, the seismic intensity distribution in the abnormal seismic region can also be estimated with high accuracy. Can be done.

Further, by the first computing program, the strength y _i and spatial coordinates of the observation point of the observed ground motion in magnitude M and a plurality of observation points i of past earthquakes (u _i, v _i), sequential new spatial coordinates (u _i, v _i) of the intensity y _i and observation point i of ground motion observed in magnitude M and observation point i during an earthquake the regression coefficient beta ₀ in _addition, by updating the beta, sequentially The accuracy of the prediction formula can be improved.

1 Receiving device (reception means)
2 PC
3 Storage device (storage means)
4 Monitor (display means)

Claims (1)

From the magnitude of past earthquakes, the strength of seismic motion observed at multiple observation points, and the spatial coordinates of the observation points, the seismic motion is obtained by obtaining the regression coefficient in the regression model using geographically weighted regression analysis in the distance attenuation formula. It is a method of constructing a strength prediction formula and estimating the distribution of the seismic intensity using the prediction formula when a new earthquake occurs .
The above prediction formula is "log Y _i = aM-log _{X-β (u i} , v _i ) X + β ₀ (u _i , v _i ) + ε _i ", where Y _i is the strength of the seismic motion at the above observation point. , M is the magnitude, X is the distance between the spatial coordinates of the observation point and the epicenter position in the past earthquakes, a, β (u _i , v _i ) and β ₀ (u _i , v _i ). Is a coefficient, ε _i is a regression error, and (u _i , v _i ) are the spatial coordinates of the above observation points.
The method, by substituting the step of calculating the W _ij by setting the value of b in formula (c) below, a W _ij calculated in the following formula (c) in formula (b) below W ( The _{process of calculating u i} , v _{i ), the process of substituting W (u i} , v _i ) calculated in the following formula (b ) into the following formula (a) to calculate β, and the following formula the beta calculated in (a) comprises a step of substituting the β (u _i, v _i) of the prediction equation,
However, in the following equation (a), X is an explanatory variable, Y is a matrix of observed values, T is a transposed matrix, and W (u _i , v _i ) is a spatial weight matrix. W (u _i , v _i ) is a spatial weighting matrix. In the following equation (c), _Wij is an element of the spatial weighting matrix, _dij is the distance between the above multiple observation points, and b is the bandwidth. A method for estimating the seismic intensity distribution.

Images