JP2014219315A - Aftershock prediction apparatus, aftershock prediction method, and program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide an aftershock prediction apparatus, an aftershock prediction method, and a program capable of estimating an aftershock probability at real time even with limited earthquake observation data right after a main shock.SOLUTION: According to the present invention, earthquake observation data by arbitrary elapsed time since at least a main shock is stored. A function based on product between a detection rate function at a detection rate of 50% at time t and including at least a parameter μ(t) indicating a magnitude and the Gutenberg-Richter law is sequentially fitted to the earthquake observation data in each learning period. Values of parameters including the parameter μ(t) varying with time t are estimated. Using an earthquake occurrence rate model into which the estimated parameter values are introduced, the temporal variation of the parameter μ(t) is estimated non-parametrically at a time of predicting an aftershock occurrence rate based on the state space law.

Description

  The present invention relates to an aftershock prediction apparatus, an aftershock prediction method, and a program.

  Conventionally, mathematical models that express the occurrence frequency of earthquakes have been proposed.

  For example, the Gudenberg-Richter rule has been discovered, which indicates that there is a proportional relationship between the magnitude of the magnitude and the logarithm of the occurrence frequency of each magnitude earthquake (aftershock).

  In 1894, Fumichi Omori published a formula (Nori Omori) indicating that the number of aftershocks after the main shock decays according to a hyperbolic function with respect to time, and in 1957 Tokuharu Utsu Announced Omori-Utsunori (Improved Omori formula), which shows that it is improved and generally attenuates according to the power law.

  In 1989, Riesenberg and Jones announced a method for predicting the probability of aftershocks at any time in the future, combining the Gutenberg-Richter rule and the Omori-Utsu rule.

  Each parameter of these earthquake occurrence models can be obtained by fitting the earthquake occurrence model to seismic observation data. Especially within one day from the main shock, individual aftershocks are well distinguished by the effect of multiple seismic waves overlapping. The resulting data is very incomplete because the aftershock detection rate is reduced. Therefore, it is necessary to consider the detection rate when fitting.

  Non-Patent Document 1 discloses that aftershock probability is estimated taking into account imperfections of seismic observation data by introducing a function of detection rate that varies parametrically according to the elapsed time from the mainshock.

Yoshihiko Ogata, Koichi Katsura “Immediate and updated forcasting of aftershock hazard”, Geophysical Research Letters, Volume 33, Issue 6

  However, the conventional method described in Non-Patent Document 1 or the like has a problem that if the earthquake observation data immediately after the mainshock is used, the detection rate function cannot be fitted and the parameters may diverge. It was.

  The present invention has been made in view of the above problems, and an aftershock prediction apparatus, an aftershock prediction method, and a program capable of estimating the aftershock probability in real time even with limited earthquake observation data immediately after the mainshock The purpose is to provide.

  In order to achieve such an object, the aftershock prediction apparatus of the present invention is an aftershock prediction apparatus including at least a storage unit and a control unit, and the storage unit stores earthquake observation data from at least the main shock to an arbitrary elapsed time. The control unit stores a function based on a product of a detection rate function including at least a parameter μ (t) indicating a magnitude of a detection rate of 50% at time t and Gutenberg-Richter rule in each learning period. A parameter estimation unit that sequentially fits observation data and estimates a parameter value including a parameter μ (t) that varies with time t, and an earthquake occurrence rate model that introduces the parameter value estimated by the parameter estimation unit Aftershock occurrence prediction means for predicting the aftershock occurrence rate using The estimation means estimates the time variation of the μ (t) non-parametrically based on a state space method.

Moreover, the aftershock prediction apparatus of the present invention is characterized in that, in the aftershock prediction apparatus described above, the time fluctuation of the μ (t) is expressed by the following discrete expression.

Moreover, the aftershock prediction apparatus of the present invention introduces the following expression as a prior distribution for the μ (t) in the aftershock prediction apparatus described above, and imposes a penalty on the time fluctuation of the μ (t). Features.
(Here, V indicates the degree of fluctuation.)

  The aftershock prediction method of the present invention is an aftershock prediction method executed in a computer having at least a storage unit and a control unit, and the storage unit stores at least earthquake observation data from the mainshock to an arbitrary elapsed time. The function of the product of the detection rate function including at least the parameter μ (t) indicating the magnitude of the detection rate of 50% at time t and the Gutenberg-Richter rule, which is executed in the control unit, is the earthquake in each learning period. A parameter estimation step for sequentially fitting observation data and estimating a parameter value including a parameter μ (t) that fluctuates according to time t, and an earthquake occurrence rate model introducing the parameter value estimated by the parameter estimation step Aftershock occurrence prediction step for predicting the aftershock occurrence rate using In the measuring apparatus, the parameter estimation step estimates the time variation of the μ (t) non-parametrically based on a state space method.

  The program of the present invention is a program for causing a computer having at least a storage unit and a control unit to store earthquake observation data from the mainshock to an arbitrary elapsed time, and for executing the time t A function based on the product of the detection rate function at least including the parameter μ (t) indicating the magnitude of the detection rate of 50% and the Gutenberg-Richter rule is sequentially fitted to the seismic observation data in each learning period, and at time t The aftershock occurrence rate is predicted using a parameter estimation step that estimates the value of the parameter including the parameter μ (t) that varies in response, and an earthquake occurrence rate model that introduces the parameter value estimated by the parameter estimation step. Aftershock occurrence rate prediction step and a program for executing In the parameter estimation step, the time variation of the μ (t) is estimated non-parametrically based on the state space method.

  The present invention also relates to a recording medium, wherein the program described above is recorded.

  According to the present invention, there is an effect that the aftershock probability can be estimated in real time even with limited earthquake observation data immediately after the mainshock. More specifically, the parametric function of the detection rate proposed in the prior art was assumed to be a parametric function in which the detection rate of earthquakes gradually increased, so that earthquake observations within one day after the occurrence of a large mainshock occurred When there is fluctuation in the detection rate in the data, there is a problem that the fitting may not be performed flexibly and the parameter may diverge. However, according to the present invention, the time variation of the detection rate is reduced by non-parametric. Therefore, even inhomogeneous seismic observation data immediately after the main shock can be optimally fitted.

  In addition, according to the present invention, the detection rate function with time as a variable is discretized and estimated by the state space method, so that adaptive estimation can be performed on data without assuming a specific function form. There is an effect that becomes possible.

FIG. 1 is a graph showing the Gutenberg-Richter rule. FIG. 2 is a graph showing a detection rate function based on a cumulative distribution function of a normal distribution. FIG. 3 is a block diagram illustrating an example of an aftershock prediction apparatus 100 to which the present exemplary embodiment is applied. FIG. 4 is a flowchart showing an example of processing of the aftershock prediction apparatus 100 in the present embodiment. FIG. 5 is a diagram showing the epicenter distribution before and after the Tohoku-Pacific Ocean Earthquake. FIG. 6 is a diagram showing magnitude time-series data before and after the Tohoku-Pacific Ocean Earthquake. FIG. 7 is a graph obtained by fitting Equation (6) based on the product of the detection rate function and the Gutenberg-Richter rule to earthquake observation data. FIG. 8 is a graph showing the value of μ (t) obtained in each learning period (3 hours, 6 hours, 12 hours, and 24 hours from the main shock). FIG. 9 is a diagram showing the relationship between the normalized histogram of magnitudes observed in each period (3-6h, 6-12h, 12-24h) and the estimated instantaneous magnitude distribution at the time of each earthquake. is there. FIG. 10 is a diagram showing a result of predicting the frequency of aftershocks having a magnitude of 5 or more in the period from 1 day to 30 days using the parameters estimated in the learning period until 1 day after the main shock. FIG. 11 is a diagram showing the results of predicting the frequency of aftershocks with a magnitude of 6 or more in the period from 1 day to 30 days, using the parameters estimated in the learning period until 1 day after the main shock. FIG. 12 is a diagram showing a result of predicting the frequency of aftershocks having a magnitude of 5 or more using the parameters estimated in the short learning period. FIG. 13 is a diagram illustrating a result of predicting the frequency of aftershocks having a magnitude of 6 or more in the subsequent period using parameters estimated in a short learning period.

  Embodiments of an aftershock prediction apparatus, an aftershock prediction method, and a program according to the present invention will be described below in detail with reference to the drawings. Note that the present invention is not limited to the embodiments.

  Hereinafter, the outline of the embodiment of the present invention will be described, and then the details of the configuration and processing of the present embodiment will be described.

[Outline of Embodiment of the Present Invention]
As an outline of the embodiment of the present invention, the background of the development of the present embodiment, prior art examples, and the outline of the embodiment of the present invention will be described in order below.

[background]
Numerous aftershocks occur after a major earthquake. In particular, during the day after the main shock, there is a high probability of a large aftershock that can cause additional damage, and about half of all large aftershocks occur during this period. Therefore, it is desirable to predict aftershocks as soon as possible after the mainshock.

  On the other hand, it is very difficult to make an accurate prediction immediately after the mainshock. In fact, the Japan Meteorological Agency makes a prediction one day after the mainshock. This is because it is difficult to determine the location, time, and magnitude of the occurrence of each earthquake due to the overlap of the main shock and many aftershocks, and only a small portion of the actual earthquakes can be observed. is there.

  In order to overcome this difficulty, it has been proposed to estimate and predict aftershock activity that actually occurred, including aftershocks from observations, taking into account the detection rate of aftershocks (Non-patent Document 1).

  In this prior art example (Non-Patent Document 1), Miyagi in 2003 was able to estimate the aftershock detection rate first and combine it with an existing prediction model to make highly accurate predictions immediately after the mainshock. It is shown as an example of aftershock activity of the Ken-Oki earthquake. Below, the prior art example of a nonpatent literature 1 is demonstrated, and the improvement method by this Embodiment is demonstrated after that.

[Prior art example]
In the prior art example (Non-Patent Document 1), a model of Reasenberg-Jones (Science 1989, Vol. 243, no. 4895, pp. 1173-1176) is used as a model of an aftershock that actually occurred. In this model, the aftershock rate of magnitude M at the elapsed time t from the mainshock is given by the following equation.
Here, K, p, c are parameters of Omori-Utsu rule (improved Omori formula) representing aftershock attenuation, and b is a Gutenberg-Richter rule parameter representing magnitude frequency. Here, FIG. 1 is a graph showing the Gutenberg-Richter rule. The horizontal axis represents magnitude (M), and the vertical axis represents the logarithm of occurrence frequency n (M).

  As shown in FIG. 1, the Gutenberg-Richter rule indicates that there is a proportional relationship between the magnitude of magnitude (M) and the logarithm of the occurrence frequency of each magnitude (n (M)). M) = a−bM (see the following formula (5)). That is, the slope of the graph of FIG. 1 is −b, and the intercept of the vertical axis is a.

From the model of equation (1), the occurrence rate of aftershocks of magnitude _Mf or higher is given by the following equation.

In order to estimate the parameters of the actual aftershock model (1) from the observed aftershock data, an earthquake detection rate is introduced. The aftershock detection rate depends on the magnitude, such that the smaller the magnitude is, the lower the magnitude is, and the higher the magnitude is, and the relationship is modeled by the following cumulative distribution function of normal distribution.
Here, the parameter μ represents the magnitude with a detection rate of 50%, and σ represents the width of the magnitude at which the earthquake is partially observed. Here, FIG. 2 is a graph showing a detection rate function based on a cumulative distribution function of a normal distribution. The horizontal axis indicates the magnitude, and the vertical axis indicates the detection rate.

  As shown in FIG. 2, μ is a magnitude (in this example, about μ = 4) when the detection rate is 0.5 (50%).

  In addition to the dependency on magnitude, the aftershock detection rate also has a dependency on time. In other words, the detection rate is very low after the mainshock, and has a time dependence that returns to the normal value over time. To express it, let the parameter μ be a function of time μ (t).

In the prior art example by Ogata and Katsura (Non-Patent Document 1), a model of the parametric time-dependent detection rate function μ (t) with parameters of μ ₀ , μ ₁ , α, γ is proposed as follows: Has been.

The parameters of the detection rate function can be estimated by fitting from the magnitude data of the observed aftershocks. In model (1), it is assumed that the distribution of magnitudes of actual earthquakes follows the Gutenberg-Richter rule:

The probability distribution of the magnitude of the detected earthquake is proportional to the product of the magnitude distribution (5) and the detection rate (3) of the earthquake that actually occurred, and is given by the following equation.

  With the above procedure, the detection rate that varies with time can be estimated using the detection rate function.

Next, the remaining parameters K, p, and c of the prediction model (1) are estimated using the estimated detection rate. The occurrence rate of the detected earthquake is expressed by the following equation which is the product of the probability intensity (1) of an actually occurring earthquake and the detection rate function (3).

From this, the logarithmic likelihood for K, p, c of Omori-Utsunori is expressed by the following equation, and the parameter value is obtained in the prior art example by maximizing this.

[Method according to this embodiment]
In the above-described prior art example (Non-Patent Document 1), a parametric statistical model is used as shown in Equation (4) in order to estimate the detection rate of aftershocks that fluctuate over time. However, when trying to estimate the parameters of this model from a short aftershock time series, especially when only the seismic observation data for a limited time immediately after the main shock is available, the parameters diverge and the optimal parameter values are obtained. I couldn't.

  As a result of intensive studies, the present inventors have found that it is appropriate to introduce a detection rate estimation method based on the state space method without using a parametric statistical model in order to estimate a stable detection rate. I came to recall.

  In the present embodiment, μ (t) is estimated using the state space method in order to stably estimate the detection rate function. By using the state space method, adaptive estimation can be performed on seismic observation data without assuming a specific function form.

In this embodiment, in order to simplify the analysis, μ (t), which is a time continuous function, is discretized as the following equation. Here, t _i is the time of occurrence of the detected aftershock to i-th.

Thus, the estimation problem of μ (t) is a finite number of parameters μ = [μ ₁ , μ ₂ ,. . . , Μ _N ] ^T becomes a problem to estimate.

In the state space method, estimation is performed by imposing a penalty on the variation of μ (t). For this purpose, a prior distribution for μ (t) is introduced as follows:

Here, the initial distribution is P (μ ₁ , μ ₂ ) = const (constant). As a result, the variation of μ (t) is limited, and a smooth function optimal for prediction can be estimated. The degree of variation is determined by the parameter V. Then, from the given μ, the magnitude sequence M = [M ₁ , M ₂ ,. . . , M _N ] ^T occurrence probability is given by the following equation from equation (6).

At this time, the target of the estimation is to maximize the posterior probability P _{β, σ, V} (μ | M) of the parameter μ from the magnitude sequence M under the optimal estimation of the parameters V, β, σ ^. Is to find. The posterior probability is given by the following equation from the Bayes formula.

To maximize the posterior distribution, μ ^* is found by the Newton method or the like. However, these estimations greatly depend on the estimation of the parameters V, β, and σ. The posterior probabilities P (V, β, σ | M) of these parameters are also given by the following equation from the Bayes formula.

  The prior distribution P (V, β, σ) may be determined from past data. As an example, it is assumed that P (V, β, σ) is a normal distribution of β. V, β, and σ that maximize P (V, β, σ | M) are obtained by fitting using the EM method (expected value maximization method) or the like.

  The above is the outline of the present embodiment. In the above, one embodiment in which the present invention is embodied is described as an outline of the present embodiment, and the present invention is not limited to the embodiment described above or described later. In the following embodiment, in order to realize the above-described embodiment by cooperation of hardware and software, details of the configuration and processing of the embodiment will be described.

[Configuration of aftershock prediction device]
Next, the structure of the aftershock prediction apparatus of this Embodiment is demonstrated with reference to FIG. FIG. 3 is a block diagram showing an example of an aftershock prediction apparatus 100 to which the present embodiment is applied, and conceptually shows a part related to the present embodiment in the configuration.

  As shown in FIG. 3, the aftershock prediction apparatus 100 according to the present embodiment generally includes at least a control unit 102 and a storage unit 106, and in the present embodiment, further includes an input / output control interface unit 108 and communication control. An interface unit 104 is provided. Here, the control unit 102 is a CPU or the like that comprehensively controls the entire aftershock prediction apparatus 100. The communication control interface unit 104 is an interface connected to a communication device (not shown) such as a router connected to a communication line or the like, and the input / output control interface unit 108 is connected to the input device 112 or the output device 114. The interface to be connected. The storage unit 106 is a device that stores various databases and tables. Each part of these aftershock prediction apparatus 100 is connected so that communication is possible via arbitrary communication paths. Further, the aftershock prediction apparatus 100 is communicably connected to the network 300 via a communication device such as a router and a wired or wireless communication line such as a dedicated line.

  Various databases and tables (such as the real-time observation file 106a, the parameter file 106b, and the prediction result file 106c) stored in the storage unit 106 are storage means such as a fixed disk device. For example, the storage unit 106 stores various programs, tables, files, databases, web pages, and the like used for various processes.

Among these components of the storage unit 106, the real-time observation file 106a is an earthquake observation data storage unit that stores earthquake observation data. The real-time observation file 106a stores at least data from the mainshock to an arbitrary elapsed time, and preferably the earthquake observation data is updated in real time according to the elapsed time from the mainshock. For example, the control unit 102 of the aftershock prediction device 100 downloads earthquake observation data (JMA catalog data, JMA centralized epicenter data, etc.) updated from the external system 200 such as the Japan Meteorological Agency via the network 300, The earthquake observation data stored in the real-time observation file 106a may be updated. In addition, the control unit 102 of the aftershock prediction device 100 acquires seismic waveform data in real time from the input device 112 such as a seismometer or from an external system 200 such as a seismometer constituting the seismic observation network, and the earthquake occurrence time is obtained. The earthquake observation data may be acquired by determining the epicenter and calculating the magnitude, and updating the earthquake observation data in the real-time observation file 106a.

The parameter file 106b is parameter storage means for storing various parameters (coefficient values and the like). For example, the parameter file 106b includes Gutenberg-Richter rule parameters (for example, β (= bln10) in the equation (5)), discretized μ (t) parameters (for example, the index i in the equation (9)) Association data of elapsed time t, μ _i (i = 1, 2,..., N)), parameters relating to aftershock attenuation (for example, parameters K, p, c relating to Omori-Utsu rule in equation (1)) , Parameters of the detection rate function relating to magnitude (cumulative distribution function of normal distribution, etc.) (for example, σ in equation (3)), parameters relating to Gaussian distribution (for example, P (μ ₁ , μ ₂ ) in equation (10), V ) And the like (initial values, update values, etc.) may be stored.

  The prediction result file 106c is a prediction result storage unit that stores the prediction result of the aftershock occurrence rate by the aftershock occurrence rate prediction unit 102b. For example, the prediction result file 106c may include an aftershock with a magnitude of 5 or more within 3 days based on a calculation result (for example, a model of the earthquake occurrence rate (Omori / Utsunori), Reasenberg-Jones model, etc.) in which parameter values are introduced. The probability of occurrence (such as XX%) may be stored.

  The input / output control interface unit 108 controls the input device 112 and the output device 114. Here, as the output device 114, in addition to a monitor (including a home TV), a speaker can be used (hereinafter, the output device 114 may be described as a monitor). In addition to the keyboard, mouse, microphone, and the like, the input device 112 may be a means for measuring earthquake motion such as a seismometer or a three-axis acceleration sensor. A plurality of input devices 112 may be connected to the aftershock prediction device 100.

  In FIG. 3, the control unit 102 includes an internal memory for storing a control program such as an OS (Operating System), a program defining various processing procedures, and necessary data. And the control part 102 performs the information processing for performing various processes by these programs. The control unit 102 includes a parameter estimation unit 102a and an aftershock occurrence rate prediction unit 102b in terms of functional concept. Schematically, the parameter estimation unit 102a estimates a parameter value by applying a mathematical model to an actual measurement value (earthquake observation data up to the present), whereas the aftershock occurrence rate prediction unit 102b estimates an estimated parameter value. It has a function to predict the rate of future aftershocks using a mathematical model that introduces.

  That is, the parameter estimation unit 102a sequentially fits a function based on the product of the detection rate function including at least the parameter μ (t) and the Gutenberg-Richter rule to the seismic observation data in each learning period, and varies according to the time t. Parameter estimation means for estimating the value of the parameter including the parameter μ (t). The parameter estimation unit 102a stores the estimated parameter value and the parameter value updated based on the newly acquired earthquake observation data in the parameter file 106b.

Here, in the present embodiment, the parameter estimation unit 102a estimates the time variation of μ (t) non-parametrically based on the state space method. More specifically, the parameter estimation unit 102a may discretize a function of μ (t) having time as a variable based on the state space method using the following expression. For example, the parameter estimation unit 102a may use a finite number (N) of magnitude series M = [M ₁ , M ₂ ,. . . , M _N ] ^T , the corresponding finite number (N) of parameters μ = [μ ₁ , μ ₂ ,. . . , Μ _N ] ^T.
Here, i is the index of the detected earthquake, and t is the elapsed time from the mainshock.

Here, the parameter estimation unit 102a may introduce the following expression as a prior distribution for μ (t). As a result, it is possible to estimate by imposing an appropriate penalty for the time variation of μ.
(Here, V indicates the degree of fluctuation.)

  Here, the parameter estimation unit 102a may perform fitting in two stages in order to estimate various parameters. That is, in the first stage, the parameter estimation unit 102a calculates a function that depends on the magnitude (formula (6) or the like representing the magnitude distribution based on the product of the detection rate function including at least the parameter μ (t) and the Gutenberg-Richter rule). The parameters of the function depending on the magnitude may be estimated for each learning period by fitting to the seismic observation data in each learning period. In this first stage, μ (t) is treated as a constant.

  Then, in the second stage, the parameter estimation unit 102a has obtained a time-dependent function (such as a function relating to the decay of the aftershock occurrence probability according to the time-dependent detection rate) from the occurrence of the mainshock. By fitting the seismic observation data up to the time, a parameter related to the time-dependent function may be estimated. For example, the parameter estimation unit 102a calculates V, β, σ that maximizes the posterior distribution P (V, β, σ | M) from the prior distribution P (V, β, σ) determined from past earthquake observation data. Alternatively, it may be obtained by fitting using the EM method (expected value maximization method) or the like.

  The aftershock occurrence rate predicting unit 102b is an aftershock occurrence rate predicting unit that predicts the aftershock occurrence rate using an earthquake occurrence rate model in which the parameter values estimated by the parameter estimating unit 102a are introduced. For example, the aftershock occurrence predicting unit 102b graphs an earthquake occurrence rate model (Omori / Utsunori, Reaseenberg-Jones model, etc.) that introduces the estimated parameter values, and graphs the future aftershock occurrence rate. Alternatively, an aftershock occurrence rate of magnitude 5 or more may be calculated by integration or the like within 3 days. Note that the aftershock occurrence rate prediction unit 102b may store the prediction result data in the prediction result file 106c, or may output (display or the like) the output result to an output device 114 such as a monitor. Etc.).

  The above is an example of the structure of the aftershock prediction apparatus 100 in this Embodiment. Note that the aftershock prediction apparatus 100 may be connected to the external system 200 via the network 300. In this case, the communication control interface unit 104 performs communication control between the aftershock prediction device 100 and the network 300 (or a communication device such as a router). That is, the communication control interface unit 104 has a function of communicating data with other terminals via a communication line. The network 300 has a function of connecting the aftershock prediction apparatus 100 and the external system 200 to each other, and is, for example, the Internet.

  The external system 200 is connected to the aftershock prediction apparatus 100 via the network 300 and provides an external database regarding various parameters, a program for causing the connected information processing apparatus to execute the aftershock prediction method, and the like. It has a function.

  Here, the external system 200 may be configured as a WEB server, an ASP server, or the like. Further, the hardware configuration of the external system 200 may be configured by an information processing apparatus such as a commercially available workstation or a personal computer and its attached devices. Each function of the external system 200 is realized by a CPU, a disk device, a memory device, an input device, an output device, a communication control device, and the like in the hardware configuration of the external system 200 and a program for controlling them.

  This is the end of the description of the configuration of the present embodiment.

[Process of aftershock prediction apparatus 100]
Next, an example of processing of the aftershock prediction apparatus 100 according to the present embodiment configured as described above will be described in detail below with reference to FIG. FIG. 4 is a flowchart showing an example of processing of the aftershock prediction apparatus 100 in the present embodiment.

As shown in FIG. 4, first, the parameter estimation unit 102a acquires the latest earthquake observation data (step SA-1). For example, the parameter estimation unit 102a downloads, from the external system 200 such as the Japan Meteorological Agency, the seismic observation data (JMA catalog data, Meteorological Agency unified seismic source data, etc.) that is updated as needed from the external system 200 such as the Japan Meteorological Agency. May store the latest earthquake observation data. In addition, the parameter estimation unit 102a acquires seismic waveform data in real time from the input device 112 such as a seismometer or the external system 200 that constitutes the seismic observation network, and determines the time and location of the earthquake and calculates the magnitude. The earthquake observation data may be acquired by specifying the earthquake by performing.

Then, the parameter estimation unit 102a in FIG. 3 performs a process of fitting the mathematical model to the earthquake observation data (step SA-2). For example, the parameter estimation unit 102a calculates the following equation (6) based on the product of the detection rate function including at least the parameter μ (t) and the Gutenberg-Richter rule, the magnitude sequence M = [ M ₁ , M ₂ ,. . . , M _N ] ^T to fit each parameter μ = [μ ₁ , μ ₂ ,. . . , Μ _N ] ^T.

At this time, the parameter estimation unit 102a limits the variation of μ (t) to follow a Gaussian distribution. That is, the parameter estimation unit 102a performs estimation while imposing a penalty on variation in μ (t). The prior distribution for μ (t) is expressed by the following equation.
The initial distribution P (μ ₁ , μ ₂ ) is constant, and V indicates the degree of variation.

The magnitude sequence M = [M ₁ , M ₂ ,. . . , M _N ] ^T occurrence probability is expressed by the following equation.

The parameter estimation unit 102a optimizes the estimation of the parameters V, β, _{and σ,} and obtains μ = μ ^* that maximizes the posterior probability P _{β, σ, V} (μ | M) of the parameter μ from the magnitude sequence M. Alternatively, the Bayesian method may be used. The posterior probability is expressed by the following equation from the Bayes formula.

The parameter estimation unit 102a may obtain μ ^* by the Newton method or the like in order to maximize the posterior distribution. On the other hand, in order to optimize the estimation of the parameters V, β, and σ, the parameter estimating unit 102a maximizes the posterior probability P (V, β, σ | M) of the following equation based on the Bayes method. Is obtained by fitting using the EM method (expected value maximization method) or the like. Here, P (V, β, σ) is assumed to be a normal distribution of β.

  Returning to FIG. 4 again, the parameter estimation unit 102a in FIG. 3 updates the parameter value newly obtained in step SA-2 by storing it in the parameter file 106b (step SA-3).

  Then, the aftershock occurrence rate prediction unit 102b substitutes the parameter value updated from the parameter file 106b for the earthquake occurrence rate model, and predicts the aftershock occurrence rate based on the earthquake occurrence rate model (step SA-4). For example, the aftershock occurrence prediction unit 102b graphs an earthquake occurrence rate model (Omori / Utsunori, Reasenberg-Jones model, etc.) that introduces estimated parameter values, and aftershocks from the main shock to the future The occurrence rate may be represented on a graph. Further, the aftershock occurrence rate prediction unit 102b may calculate the probability that an aftershock of magnitude 5 or more will occur within 3 days by integration or the like from the earthquake occurrence rate model in which the estimated parameter values are introduced. Note that the aftershock occurrence rate prediction unit 102b may store the prediction result data in the prediction result file 106c, or may output (display or the like) the output result to an output device 114 such as a monitor. Etc.).

  And the control part 102 of the aftershock prediction apparatus 100 determines whether a process is complete | finished (step SA-5). For example, the control unit 102 of the aftershock prediction apparatus 100 may determine that the process is finished when a certain period of time has passed since the mainshock, and the process is performed when the end of the process is instructed by the user via the input unit 112 May be terminated.

  When it determines with not complete | finishing a process (step SA-5, No), the control part 102 of the aftershock prediction apparatus 100 repeats the process of step SA-1-SA-5 mentioned above. On the other hand, when it determines with complete | finishing a process (step SA-5, Yes), the control part 102 of the aftershock prediction apparatus 100 complete | finishes a process.

  The above is an example of the process of the aftershock prediction apparatus 100 in this embodiment.

[Example]
Next, an embodiment in which the above-described embodiment is applied to aftershocks data of the 2011 Tohoku-Pacific Ocean Earthquake (M9.0) will be described below with reference to FIGS.

  Here, FIG. 5 is a diagram showing the epicenter distribution before and after the Tohoku-Pacific Ocean Earthquake. The epicenter of the mainshock is represented by an asterisk. FIG. 6 is a diagram showing magnitude time-series data before and after the Tohoku-Pacific Ocean Earthquake. The horizontal axis is the time when the occurrence time of the mainshock is 0, and the vertical axis indicates the magnitude. These seismic observation data were obtained from PDE / NEIC catalog data.

  Fig. 7 is a graph of fitting the equation (6) based on the product of the detection rate function and the Gutenberg-Richter rule to μ (t) constant for the earthquake data before the Tohoku-Pacific Ocean Earthquake. is there. The horizontal axis is magnitude, and the vertical axis is frequency of occurrence. The solid line represents an expression based on the product of the detection rate function and the Gutenberg-Richter rule, and the black circle represents the magnitude frequency. Since the model and the observation are in good agreement, the validity of the detection rate function (3) is known.

  FIG. 8 is a graph showing the value of μ (t) estimated by the state space method for earthquake observation data in each learning period (3 hours, 6 hours, 12 hours, and 24 hours from the main shock). Black spots represent each observed earthquake. The dotted line in FIG. 8 shows the transition of the value of μ (t) obtained in the learning period of 3 hours, the alternate long and short dashed line is 6 hours, the alternate long and two short dashes line is 12 hours, and the solid line is 24 hours. Gray circles represent each detected earthquake. As shown in FIG. 8, it was found that the detection rate can be estimated appropriately even when the fluctuation of the detection rate immediately after the mainshock fluctuates.

  Here, FIG. 9 shows the relationship between the normalized histogram of the observed magnitude and the estimated instantaneous magnitude distribution at the time when each earthquake occurred in each period (3-6h, 6-12h, 12-24h). FIG. The horizontal axis represents magnitude, and the vertical axis represents probability density. The black circle shows the normalized histogram of the magnitude of each observed earthquake, the gray curve shows the magnitude distribution at the moment when each earthquake occurred, and the black curve is the average curve. In this way, it was confirmed that the estimated value at the moment when each earthquake occurred was in good agreement with the measured value.

  Next, using the estimated detection rate, parameters of the prediction model (Equation (1)) were estimated and a prediction experiment was performed. FIG. 10 and FIG. 11 predicted the frequency of aftershocks of magnitude 5 or more and magnitude 6 or more, respectively, from 1 day to 30 days, using the parameters estimated in the learning period until 1 day after the main shock. It is a figure which shows a result. The solid line is the function based on the product of the detection rate function estimated by this embodiment and the Gutenberg-Richter rule. The curve obtained by fitting Omori and Utsunori directly without considering the detection rate is shown.

  FIG. 12 and FIG. 13 are diagrams showing the results of predicting the subsequent period, the magnitude of magnitude 5 or more, and the frequency of aftershocks of 6 or more, respectively, using parameters estimated in a shorter learning period. The one-dot line indicates the prediction result obtained in the learning period of 3 hours, the broken line indicates 6 hours, the two-dot chain line indicates 12 hours, and the solid line indicates the learning period of 24 hours. As a result, it was confirmed that the prediction according to the present embodiment and the actual observation were in good agreement.

[Other embodiments]
Although the embodiments of the present invention have been described so far, the present invention is not limited to the above-described embodiments, but can be applied to various different embodiments within the scope of the technical idea described in the claims. It may be implemented.

  For example, although the case where the aftershock prediction apparatus 100 performs processing in a stand-alone form has been described as an example, the aftershock prediction apparatus 100 performs processing in response to a request from the client terminal, and returns the processing result to the client terminal. You may make it do.

  In addition, among the processes described in the embodiment, all or part of the processes described as being automatically performed can be performed manually, or the processes described as being performed manually can be performed. All or a part can be automatically performed by a known method.

  In addition, unless otherwise specified, the processing procedures, control procedures, specific names, information including registration data for each processing, parameters such as search conditions, screen examples, and database configurations shown in the above documents and drawings Can be changed arbitrarily.

  In addition, regarding the aftershock prediction apparatus 100, each illustrated component is functionally conceptual and does not necessarily need to be physically configured as illustrated.

  For example, the processing functions provided in each device of the aftershock prediction device 100, in particular, the processing functions performed by the control unit 102, all or any part thereof are interpreted and executed by a CPU (Central Processing Unit) and the CPU. It may be realized by a program to be executed, or may be realized as hardware by wired logic. The program is recorded on a recording medium to be described later, and is mechanically read by the aftershock prediction apparatus 100 as necessary. That is, in the storage unit 106 such as a ROM or an HD, a computer program for performing various processes by giving instructions to the CPU in cooperation with an OS (Operating System) is recorded. This computer program is executed by being loaded into the RAM, and constitutes a control unit in cooperation with the CPU.

  The computer program may be stored in an application program server connected to the aftershock prediction apparatus 100 via an arbitrary network 300, and may be downloaded in whole or in part as necessary. It is.

  In addition, the program according to the present invention may be stored in a computer-readable recording medium, and may be configured as a program product. Here, the “recording medium” means a memory card, USB memory, SD card, flexible disk, magneto-optical disk, ROM, EPROM, EEPROM, CD-ROM, MO, DVD, and Blu-ray (registered trademark). It includes any “portable physical medium” such as Disc.

  The “program” is a data processing method described in an arbitrary language or description method, and may be in any format such as source code or binary code. The “program” is not necessarily limited to a single configuration, but is distributed in the form of a plurality of modules and libraries, or in cooperation with a separate program represented by an OS (Operating System). Including those that achieve the function. Note that a well-known configuration and procedure can be used for a specific configuration for reading a recording medium, a reading procedure, an installation procedure after reading, and the like in each device described in the embodiment.

  Various databases and the like (real-time observation file 106a, parameter file 106b, prediction result file 106c, etc.) stored in the storage unit 106 are memory devices such as RAM and ROM, fixed disk devices such as hard disks, flexible disks, and optical disks. And stores various programs, tables, databases, web page files, and the like used for various processes and website provision.

  The aftershock prediction apparatus 100 may be configured as an information processing apparatus such as a known personal computer or workstation, or may be configured by connecting any peripheral device to the information processing apparatus. Further, the aftershock prediction apparatus 100 may be realized by installing software (including programs, data, and the like) that causes the information processing apparatus to realize the method of the present invention.

  Furthermore, the specific form of distribution / integration of the devices is not limited to that shown in the figure, and all or a part of them may be functional or physical in arbitrary units according to various additions or according to functional loads. Can be distributed and integrated. That is, the above-described embodiments may be arbitrarily combined and may be selectively implemented.

  As described above in detail, according to the present invention, an aftershock prediction device, an aftershock prediction method, and a program capable of estimating the aftershock probability in real time even with limited earthquake observation data immediately after the mainshock It is extremely useful in various fields such as disaster prevention and earthquake prediction, such as the announcement of aftershock predictions by the Japan Meteorological Agency, etc., and the service of distributing disaster information predictions by private companies for a fee.

DESCRIPTION OF SYMBOLS 100 Aftershock prediction apparatus 102 Control part 102a Parameter estimation part 102b Aftershock incidence prediction part 104 Communication control interface part 106 Storage part 106a Real time observation file 106b Parameter file 106c Prediction result file 108 Input / output control interface part 112 Input device 114 Output device 200 External System 300 network

Claims (5)

An aftershock prediction device having at least a storage unit and a control unit,
The storage unit
Store at least earthquake observation data from the mainshock to any elapsed time,
The control unit
A function based on the product of the detection rate function at least including the parameter μ (t) indicating the magnitude of the detection rate of 50% at time t and the Gutenberg-Richter rule is sequentially fitted to the seismic observation data in each learning period, and time t Parameter estimation means for estimating the value of a parameter including a parameter μ (t) that varies according to
An aftershock occurrence predicting means for predicting an aftershock occurrence rate using an earthquake occurrence rate model in which the value of the parameter estimated by the parameter estimation means is introduced;
Aftershock prediction device with
The parameter estimation means includes
An aftershock prediction apparatus characterized by estimating the time fluctuation of the μ (t) non-parametrically based on a state space method. The aftershock prediction apparatus according to claim 1,
The time variation of μ (t) is
Represented by the following equation, which is discretized based on the state space method
Aftershock prediction device characterized by. In the aftershock prediction device according to claim 2,
The following formula is introduced as a prior distribution for the μ (t), and a penalty is imposed on the time variation of the μ (t).
(Here, V indicates the degree of fluctuation.)
Aftershock prediction device characterized by. An aftershock prediction method executed in a computer having at least a storage unit and a control unit,
The storage unit
Store at least earthquake observation data from the mainshock to any elapsed time,
Executed in the control unit,
A function based on the product of the detection rate function including at least the parameter μ (t) indicating the magnitude of the detection rate of 50% at time t and the Gutenberg-Richter rule is sequentially fitted to the seismic observation data in each learning period, and at time t A parameter estimation step for estimating a value of a parameter including a parameter μ (t) that varies in response;
An aftershock occurrence rate predicting step for predicting an aftershock occurrence rate using an earthquake occurrence rate model in which the value of the parameter estimated in the parameter estimation step is introduced;
In aftershock prediction equipment including
The parameter estimation step includes
A method for predicting aftershocks, wherein the time variation of μ (t) is estimated non-parametrically based on a state space method. A program for causing a computer having at least a storage unit and a control unit to store earthquake observation data from at least a main shock to an arbitrary elapsed time,
In the control unit,
A function based on the product of the detection rate function including at least the parameter μ (t) indicating the magnitude of the detection rate of 50% at time t and the Gutenberg-Richter rule is sequentially fitted to the seismic observation data in each learning period, and at time t A parameter estimation step for estimating a value of a parameter including a parameter μ (t) that varies in response;
An aftershock occurrence rate predicting step for predicting an aftershock occurrence rate using an earthquake occurrence rate model in which the value of the parameter estimated in the parameter estimation step is introduced;
A program for executing
The parameter estimation step includes
A non-parametric estimation of the time variation of μ (t) based on a state space method.