Research Memorandum No. 916

June 14, 2004

Preliminary Analysis of Observations on the Ultra-Low
Frequency Electric Field in a Region around Beijing

Jiancang ZhuangDavid Vere-Jones
Huaping Guan
Yosihiko Ogata
and
Li Ma

Preliminary Analysis of Observations on the Ultra-Low
Frequency Electric Field in a Region around Beijing
Jiancang Zhuang1 , David Vere-Jones2 , Huanping Guan 3
Yosihiko Ogata1 and Li Ma3

Abstract — This paper presents a preliminary analysis of observations on ultra-low frequency groundelectric signals from stations operated by the China Seismological Bureau over the last 20 years. A brief
description of the instrumentation, operating procedures and data processing procedures is given. The
data analysed consists of estimates of the total strengths (cumulated amplitudes) of the electric signals
during 24 hour periods. The thresholds are set low enough so that on most days a zero observationis returned. Non-zero observations are related to electric and magnetic storms, occasional man-madeelectrical effects, and, apparently, some pre-, co-, or post-seismic signals. The main purpose of the analysis
is to investigate the extent that the electric signals can be considered as pre-seismic in character. For
this purpose the electric signals from each of five stations are jointly analyzed with the catalogue of local
earthquakes within circular regions around the selected stations. A version of Ogata’s Lin-Lin algorithm isused to estimate and test the existence of a pre-seismic signal. This model allows the effect of the electric
signals to be tested, even after allowing for the effects of earthquake clustering. It is found that, although
the largest single effect influencing earthquake occurrence is the clustering tendency, there remains asignificant preseismic component from the electrical signals. Additional tests show that the apparent
effect is not post-seismic in character, and persists even under variations of the model and the time
periods used in the analysis. Samples of the data are presented, and the full data sets have been madeavailable on local websites.
Key words: Ultra-Low frequency electric signal, earthquake risk, Hawkes’ self-exciting and mutually
exciting model, Beijing

1

Introduction

Systematic studies on the ground electric field as a source of possible earthquake precursors were started in
China following reports of substantial electric anomalies before the MS 7.8 Tangshan earthquake in 1976.
Since that time more than 100 stations have been built in China for observing the ground electric field.
Initially, observations were made over the whole frequency band. Experience showed that observations
on specific frequency bands, particularly the frequencies from 0.1 to 20 Hz, were most effective and best
able to avoid contamination by industrial noise. In 1981, construction began on an ultra-low-frequency(ULF) observational network in the Beijing region, using a frequency band of 0.1-10 Hz. The first stations
from this network started operation in 1982, while others were added during the ensuing decade. The
main purpose of this paper is to provide a preliminary statistical analysis of accumulated data from thestations in this network.
Earlier studies of the electric signals data, of a more informal character and involving stations fromthe whole of China, were made by Chen and Xu (1994), and Guan et al. (1995, 1996, 1999, 2000). Theysuggested that anomalous fluctuations of the ultra-low frequency electric field may occur up to twenty
days before a relatively large earthquake in a neighborhood of the future focal region, and that the characteristics of the anomalies may be related to the distance of the recording station from the epicenter, andthe magnitude of the coming earthquake. This paper concentrates on data from five selected stations in the
Beijing network, and uses a statistical model to try to quantify the correlations between the occurrence of
anomalous fluctuations in the electric field, and the occurrence of local earthquakes. The earthquake data
1 Instituteof Statistical Mathematics,4-6-7 Minami Azabu, Minato-Ku, Tokyo 106-8569, Japan. Email:zhuangjc@ism.ac.jp
2 School of Mathematical and Computer Science, Victoria University of Wellington, P.O. Box 600, Wellington, New
Zealand.
3 Center for Analysis and Prediction, China Seismological Bureau, 63 Fuxing Road, Beijing 100036, China.

1

used in this study are listed at the end of the paper, and were extracted from the China National Catalogue. A sample of the electric field data used in this study is also listed at the end of the paper; the full set,
together with the earthquake data, is available from website http://www.ism.ac.jp/∼ogata/RM916/Data.
In the following sections, we first give a brief account of the observational equipment and operatingprocedures used by stations in the network (Section 2). This is followed in Section 3 by an outline ofthe procedures used to obtain the earthquake and electric signals data in the form used in the analysis.
Section 4 describes the main statistical model used, which is applied in Sections 5 and 6 to data fromthe five selected stations. The model is used to check the presence of both pre- and post-seismic effects,
and allows for the effects of earthquake clustering. Section 6 also outlines some supplementary analyses
designed to check internal consistency and other features. The final section sets out the main conclusionsof the paper.

2

Observation Network and Equipment

The electric signals observation network, from which the data for the present study was taken, is locatedin a broad region around Beijing: see Fig 1. It was started in 1981 with the aim of monitoring fluctuationsin the electro-magnetic radiations at frequencies in the 0.1–10 Hz (ULF) frequency band. The station
locations were chosen to avoid, as far as possible, man-made sources of interference from electric railways,
highways, underground metal pipes, power transformer substations, high-voltage power lines, areas of highenergy consumption, radio stations and power lines with ground connections. At the same time they wereselected to be close to known seismic belts or active faults. 8 stations are currently operating in thisnetwork, from which the 5 stations with the longest and highest quality records were selected for thepresent study. The locations of the selected stations, namely Langfang, Sanhe, Qingxian, Huailai and
Changli, are also shown in Figure 1. Langfang, Sanhe, and Qingxian stations were set up in 1981–1982,
Huailai in 1987, and Changli in 1990.
The equipment and observation systems are similar at all the above stations. They use the E-EM
system designed by the Provincial Seismological Bureau of Hebei. We provide only a brief description
here, referring to Chen and Xu (1994), Chen et al. (1998). for more detailed illustrations and technicalparameters of the operating systems.
The basic structure of the E-EM system is shown in Figure 2(a). It uses HBEMD-3 type sensors(see Chen et al, 1998) to detect the electric signals. To reduce the effects of polarization potential, theelectrodes of the sensor are made of high-quality stainless steel (Cr18Ni9C), and are cylinder shapedwith a height of 300mm and a radius of 4 mm; see Figure 2(b). In addition, 1000µFarad capacitors areconnected in series into the system at each station, to eliminate direct currents.
At each station, two pairs of electrodes are installed along perpendicular axes. The electrodes areburied from 6 meters to 12 meters deep, and some 40 meters apart. The electric signals detected by thesensor are transmitted to the pre-processor for preliminary noise reduction, and then to the recorder. Allwires used for signal transmission are screened by high quality metal nets, covered by water-proof pipes,
and buried 0.6–0.8 meters below the surface. They go directly to the laboratory from underground, and
it is checked that no power supply wires, communication cables or metal blocks are present nearby before
operation commences.
The pre-processor includes an impedance matching unit and filter. Its main functions are to prevent
any 50 Hz noise produced by industrial power supplies from coming into the amplifier, and to preventthe polarization potential from changing the working point of the amplifier.
The signals are then fed to the amplifier to activate the recorder pen. The recorders used (shown in
Figure 2(c)) are a modification of the DJ-1 recorder, originally designed for recording the medium-longperiod components (0–10 HZ) of ground movement with the DK-1 seismometer. These recorders andseismometers are widely used in Chinese seismological stations. The recording method is by automatic
continuous pen record.
In order to improve the power supply quality, a stabilization plant is used for the electricity supply
instead of the usual commercial power supply. Power from automatically recharged batteries is used
when the power from the stabilization plant is cut off. In addition, dual T bandpass filters, together withimpedance matching methods, are used to filter out any 50 Hz signals deriving from local power supplies
or other sources.
The resulting record sensitivity is 0.5 mV/mm, with a noise level of less than 0.15 mV. The frequency
response function of the system is shown in Figure 3. It confirms that the main frequency band is 1-10Hz.
2

Noise reduction is a crucial feature of the operational system. Considerable care is needed in selecting
the site of the station, and adjusting the sensitivity so that the smaller anomalous signals lie just above thenoise threshold. Other design features which help in noise reduction include the choice of the observationalfrequency band, the depths to which the electrodes and transmission cables are buried, the use of highquality electrodes to reduce the effects of polarization potential, use of a high input resistance, screeningof all transmission cables, filtering out of 50 Hz noise, and power-supply stabilization. Even despite these
efforts, the existing stations differ considerably in their ability to pick up or distinguish the anomalous
signals.

3

Data

3.1

Electric Signals Data

The data used for the present analysis is the list of signal strengths, as reported each day from each of theelectric signals stations to the Centre for Analysis and Prediction in Beijing. As already mentioned, thefive stations for the present study, chosen on the basis of the quality and completeness of their records,
are Huailai, Changli, Sanhe, Qingxian and Langfang. Because the stations were set up at different times,
and operated continuously over different periods, the observation time intervals for these stations aredifferent and are set out in detail in Table 1.
The signal strengths are determined on a daily basis according to the following protocol, which
was established at the beginning of the observation period, and is observed in the same manner at eachstation. Each day, the drum record for each pair of electrodes is examined for the presence of anomalies.
The threshold of the observing system is set low enough so that on most days the recording is close toa horizontal line (zero). Typically, the anomalous signals do not occur continuously throughout the day,but in episodes. A given episode may assume a variety of forms, but most commonly appears on the drum
record as a signal of roughly sinusoidal character with irregular amplitude and frequency (see Fig 4).
The measure adopted to quantify the daily signal strength is a rough estimate of the total cumulative
amplitude (duration times mean amplitude), computed as set out below.
First, for each episode in each of the two measurement directions (N-S and E-W), an average
amplitude is estimated from the drum record by taking half the maximum throw on the chart plus half
an approximate mean square value. The average amplitude of the episode is then converted to an average
electric field strength by multiplying by an instrument scaling factor (mV/mm) and dividing by the
distance between the electrodes.
Next, the average field strength for each episode is multiplied by the duration of the episode (inseconds, also taken from the drum record) to produce a measure A of the total strength of that particular
episode. If N episodes occur during the day, with associated total strengths (A1 , . . . , AN ), the total
cumulated daily strength is determined as
Atotal = A1 + A2 + · · · + AN .
This procedure is used for each of the two measurement directions (N-S and E-W), and yields thetwo components AN S , AEW of the strengths which are reported each day to CAP in Beijing.
A sample of the records of the daily strengths from one of the stations is shown in Table 2. Fromthese it can be seen that the signal strengths vary greatly from day to day. Zero values indicate thedays on which no anomalous activity could be detected; non-zero values vary from one to around one
thousand. The full set of daily strengths from the five stations is too lengthy to be listed in the paper,but is available from the ISM Website http://bemlar.ism.ac.jp/Data/.
For the purpose of the present analysis only, detailed variation of the daily strengths was ignored,
and the days were crudely classified as either possessing (1) or not possessing (0) anomalous signals,
according as to whether or not the the daily strengths satisfied the criterion
A = AN S + AEW ≥ 200
In this way, the signal data was transformed into a set of daily 0−1 values which were conventionallyregarded as occurring at 0.00 hours on the day following the signal readings. This series was then used as
one of the components in the later correlation studies. Table 1 lists for each station the numbers of days on
which earthquake events were recorded, the numbers of days on which signal events (1’s) were recorded,
and the ratio between these numbers. The last column is included as a possible indicator of the noisiness3

of the site, in terms of the number of electric signal days not directly associated with earthquakes. The
occurrences are displayed in Figure 5, together with the earthquake data for each station.

3.2

Recognition of interference signals

Experience has shown that low-sky or sky-to-ground lightning, power lines that leak electricity to theground near the station, and local domestic or industrial activities may influence electric field observations
made using the buried electrode method. Some such sources of noise can be distinguished easily fromthe anomalous signals described earlier by the differences in character of the recorded waveforms. As an
example, Figure 6 shows both waveforms typical of the anomalous signals, and a single large pulse typicalof an interference signal. Figure 6 shows the interference signals from low-sky lightning and power lines
leakages near the station.
When clearly distinguishable features of the latter kind are observed, they are excluded from thecalculation of the daily strengths. In all other cases, irrespective of their supposed source, fluctuations
above the noise level are assumed to form part of the anomalous signal, and are included in the calculationof the daily strengths.

3.3

Earthquake Data

For each station, a special sub-catalogue was prepared, consisting of events for which either(a) the epicentral distance from the station was less than 200 km, and the local magnitude was 4.0 or
greater, or(b) the epicentral distance from the station was less than 300 km, and the local magnitude was 5.0 or
greater.
The data for these subcatalogues was extracted from the main China catalogue prepared by CAP;
copies of this catalogue, if required for research purposes, can be obtained from the weblink http://wwwwdcds.seis.ac.cn (contact Professor Ma Li in case of difficulties, as the Website is currently in Chinese
only).
An epicentral plot of earthquakes with local magnitudes 4 and greater occurring within an approximately 800 km square region around Beijing is shown in Figure 1. Earthquakes falling into at least oneof the five sub-catalogues are shown with dark circles; other earthquakes are shown with light circles.
Time-magnitude plots of the earthquakes from the subcatalogues for each of the five selected stations areshown in Figure 5.
The five sub-catalogues are also condensed and summarized in Table 3. The table lists origin times,
epicentres and local magnitudes for each event, and indicates the stations for which the given event
fell within one or other of the two specified circles. For each station, the earthquake events satisfying
the prescribed criteria were coded as a sequence of occurrence times, ignoring magnitudes and othercoordinates. The two series for each station, one series for the electric signals and the other for theearthquakes, formed the inputs for the analyses described below.
None of the subcatalogues has been declustered prior to the analysis; rather, the clustering effectis itself modelled through the self-exciting term in the models discussed in the next section. Indeed, asecondary aim of the paper is to check the extent to which clustering can affect the apparent significanceof the electric signals terms (cf Michael, 1997). With the exception of Huailai station, for which the 300km
observation region includes the aftershock regions of both the Datong and Zhangbei events, the numberof direct aftershocks entering the subcatalogues is rather small, as can be verified from Figure 5.

3.4

Some features of the data

From the preceding tables and plots, several important features can be seen which may be useful to bearin mind for the subsequent analysis.
1. Much of the earthquake activity during the study period is associated with three major earthquake
groups. The first represents continuing low-level activity in an extended region around the epicentreof the 1976 MS =7.8 Tangshan earthquake. This activity continued into the early 1990’s. The othertwo major groups are mainly aftershocks of the two M = 6.4 Datong earthquakes of 1991/1992, andof the 1998 M = 6.2 Zhangbei earthquake. The whole region was quiet between these two events.
4

2. The clustering effect is particularly pronounced for the Huailai station, the study region for which
encompassed both the Datong and Zhangbei clusters, as well as partially extending into the regionof the Tangshan events. By contrast, the Changli station entered the study only in 1990, and lies
too far to the East to contain the smaller (4 ≤ M < 5) aftershocks from either the Datong or the
Zhangbei sequences. Its subcatalogue shows the least clustering of all five stations.
3. The electric signals are also very highly clustered, even on some occasions when there appear to
be no associated earthquakes. This feature is particularly pronounced for Langfang and Huailai
stations, both of which show substantially increased electric signals activity in the later part of therecord. These two stations, especially Langfang, were affected by urbanization of their immediate
environment during the study period. Urbanization has not occurred to the same extent at theother stations. The increased activity observed at the first two stations may therefore give someindication of the sort of interference effects to be expected from man-made sources.

4

Self-exciting and mutually exciting models

In this section we describe the model used for the major part of the analysis. The self-exciting andmutually exciting earthquake model was developed by Ogata and Utsu (see Ogata et al. 1982; Ogata,
1983; Utsu and Ogata, 1997) from the Hawkes process (Hawkes, 1971). It is most easily described throughits conditional intensity function,
λ(t)dt = E[N (dt)|Ht ],(1)
where Ht = {Observation history up to time t}. In essence, the conditional intensity function represents
the target process as a time-varying Poisson process with rate λ(t) conditioned by the past history. The
conditional intensity function of the combined self-exciting and mutually-exciting model (referred to as
the combined model in the rest of the paper) can be written as
λ(t) = µ + λS (t) + λE (t),(2)

where µ represents the constant background rate, λS is the self-exciting term, which models clustering
among the target events, and λE is the external excitation term, which models the contribution to therate from the external process. The self-exciting term is taken in the formXλS (t) =g(t − ti ),
(3)
ti <t

with the summation extended to all the events occurring before time t in the target process {t i : i =1, 2, · · · , n1 }, g(t) being a sum of Laguerre polynomialsg(t) = e−αtNSXp k tk .

(4)

k=0

Similarly, the external excitation term is written asX
λE (t) =h(t − ui ),

(5)

ui <t

with the summation taken over all the events occurring before time t in the process of precursor events(here the electric signals) {ui : i = 1, 2, · · · , n2 }, h(t) again being a sum of Laguerre polynomials,
h(t) = e

−βt

NE
X

q k tk .

(6)

k=0

The Laguerre polynomials are used here as a convenient family of orthogonal functions which inprinciple can be used to model the response functions to any required degree of accuracy; more detailed
discussions of the model are given in Vere-Jones and Ozaki (1982), Ogata et al (1982), Ma and Vere-Jones(1997), or the IASPEI manual, Utsu and Ogata (1997). The model was used in the earlier studies mainly
to investigate possible triggering effects between different kinds of seismicity.

5

In practice the number of terms included in the sum has to be balanced between concerns ofsensitivity and over-fitting, and is the main target of the model-selection procedures described later.
With limited data, as in the present situation, and little prior knowledge as to the likely form of the
response functions, the fitted response functions can be interpreted only as crude approximations to any
underlying physical processes.
Usually, the excitation effect in the fitted model reaches its maximum immediately or shortly afteran event occurs, then decays quickly with time and becomes negligeable after sufficient time.
If we drop out λE , the model becomes a self-exciting model; if we drop out λS , the model becomes
an externally excited model; if both λS and λE are neglected, the model reduces to a Poisson model. Allof these types of models will be used in our analysis.
We shall keep the notation {ti : i = 1, 2, · · · , n1 } for the main or target process, and {ui : i =1, 2, · · · , n2 } for the external exciting process or secondary process.
In our analysis, we first take the earthquake events as the target process and the anomalous electric
signals as the external process, to see whether the electric signals have explanatory power in relation tothe occurrence times of earthquakes. Then we reverse the roles of the two processes, take the anomaly
events as the main process and earthquakes as the external process, to see whether the electric signals
might be triggered by some mechanism following the occurrence of an earthquake.
Given a set of observation data, the parameters of the model can be estimated by maximizing the
likelihood function, which for a conditional intensity model has the standard form (analagous to that fora Poisson process)Z TX
λ(t; θ)du,(7)
log L(θ) =log λ(ti ; θ) −0

i:0≤ti ≤T

where ti denotes the occurrence time of the ith event, [0, T ] is the observation time interval, and Θ isthe vector of parameters (see, for example, Daley and Vere-Jones, 2002, Chapter 7).
The parameters for the combined model can be written asθ = (µ; α, p1 , p2 , · · · pNS ; β, q1 , q2 , · · · , qNE ).
Model selection, particularly the determination of the numbers of parameters NS and NE in the Laguerre
expansions, was carried out using the Akaike Information Criterion (AIC, see Akaike, 1974). The statistic
AIC = −2 max log L(θ) + 2kp
θ(8)

is computed for each of the models fitted to the data, where kp is the total number of fitted parameters. In
comparing models with different numbers of parameters, addition of the quantity 2kp roughly compensates
for the additional flexibility which the extra parameters provide. The model with the lowest AIC value
is taken as giving the optimal choice for forward prediction purposes.
Insofar as it depends on the likelihood ratio, the AIC can also be used as a rough guide to model
testing. As a rule of thumb, in testing a model with k + d parameters against a null model with just kparameters, we take a difference of 2 in AIC values as a rough estimate of significance at the 5% level.
If standard asymptotics were applied, such a difference would correspond to a significance level of 4.6%
when d = 1, 5% when d = 2, 4.6% when d=3, and 3.5% when d=5. Such figures give a rough guideto significance levels, but should be used conservatively, because of the relatively small sample sizes andother approximations.
One of the main advantages of the combined model in the present context is that it allows the effectof the electric signals terms to be examined even in the presence of clustering (modelled by self-excitation)in the earthquakes themselves.

5

Main Analysis

In the main analysis, we use the model format available within the Lin-Lin programme (Utsu and Ogata,1997). This requires all coefficients in the trend and polynomial expansions to be non-negative, thus
ensuring positivity of the conditional intensity, at the expense of some flexibility in the functional forms
available.
The combined model is fitted separately to data from the five selected stations in the Beijing region,
namely Huailai, Changli, Sanhe, Qingxian and Langfang. Because of data availability, the observation
time intervals for those stations are different, as summarized in Table 1.
6

Taking the earthquakes as the target events, we consider the following models to examine therelationship between the electric signals and the earthquakes:
1. Poisson process with polynomial trend (restricted to second order);2. self-exciting model, without external excitation;
3. externally excited model without self-excitation;
4. combined model of both self-exciting and external excitation terms.
In Section 6, to test whether the electric signal might be a post-seismic effect, we interchange theroles of the electric signals and the earthquakes.

5.1

Self-exciting versus Poisson models.

The purpose of this step is to determine the size and significance of the clustering effect among theearthquakes before using the self-exciting model as a base model to determine the contribution of theelectric signals in fitting the earthquake data. This step was accomplished by fitting the self-exciting
model to the earthquake data, and comparing it to the Poisson model with polynomial trend. The firstand second row blocks in Table 4 list the outputs from this analysis for each of the five selected stations.
From this table we can see that the clustering effect plays an important role in the process of theearthquake occurrences. For each station except Changli the reduction in AIC is very substantial, well
above what would be required to establish significance at around the 5% level. The effect is particularly
pronounced for Huailai, as we might expect from the clusters contained within the Huailai record. Thechanges in seismicity rate for Huailai were mainly caused by these clusters, so here and elsewhere we haverestricted the trend to a second degree polynomial with non-negative coefficients. The Changli station,by contrast, shows a relatively small degree of clustering, large enough for the self-exciting model to bepreferred to the Poisson model, but barely large enought to establish the significance of the clustering
effect.
We conclude that clustering should be taken into account for all five stations, and that it is aparticularly important feature for the stations closest to the source regions of the Datong and Zhangbei
clusters.

5.2

Externally excited versus Poisson models

We next compare the Poisson model to the externally excited models for the five data sets (see the third
row block in Table 4. We see that the electric signals reduce the Poisson AIC values by amounts which areless than the reductions due to the self-exciting terms for Huailai and Langfang stations, but comparable
to or greater than those reductions for Sanhe, Changli and Qingxian stations. These large differencessuggest that the external signals have considerable explanatory power. However, we shall see that thesevalues are somewhat inflated, being based in part on the fact that both earthquakes and signals are
clustered, so that to a degree the electric signals can act as a surrogate for the self-exciting terms.

5.3

Combined versus self-exciting models.

Finally, we compare the fits of the self-exciting and combined models. The final row block of Table 4 showsthat addition of the external excitation terms to the self-exciting terms contributes further reductions of
the AIC values for all five stations. All five reductions are substantial, even though smaller than thoseobtained by testing the electric signals models directly against the Poisson model. Note that the differencesin AIC values between the combined model and the Poisson model are much larger than the differences
between the externally excited model and the Poisson. Both effects show that the explanatory powerof the electric signals is considerably exaggerated unless the clustering terms are taken into account, as
forewarned by Michael (1997). The important conclusion, however, is that the electric signals retain asignificant explanatory power even after the clustering has been taken into account by the self-exciting
terms.
Additional insight can be obtained by examining the information gains per event. This quantity- here just the difference in log-likelihoods normalized by the number of events - is a measure of the
improvement in predictability in passing from the base model to the test model (see e.g. Vere-Jones7

(1988) or Harte and Vere-Jones (2004) for further discussion; the idea goes back to early papers by
Kagan (eg Kagan and Knopoff, 1977). The gains per earthquake event for each station are shown in
Table 5, together with the gains per signal event in passing from the self-exciting to the combined model.
It is interesting that the stations showing the smallest degree of clustering show the largest gains/event.
The low gains per electric signal for Huailai and Langfang suggest that these stations may be more subject
to interference from noisy signals than the other stations, although the low gains may also represent justa further nuisance effect of the clustering.

5.4

Fitted conditional intensity and impulse response functions

An example of the conditional intensity function of the combined model, and the associated impulse
response functions, is shown in Figure 7 for Sanhe station. The response functions are only crudely fitted
because of the rather small amount of earthquake data. The shapes of the impulse response functions
g(t) and h(t) suggest that the relative risk of earthquakes is high during the first few days after theoccurrence of an earthquake or an electric signal event, then decays rapidly with time. Note, however,the differences in scales for the self-exciting and electric signals terms. The latter is much more diffuse,
although the apparent rate of decay for both response functions may be exaggerated by the choice of
Laguerre functions, associated with exponential decay, as the function basis (see further discussion in
Section 6). The differences between the maximum absolute values of two functions g(t) and h(t) is anindication that, for this station at least, the immediate increase in risk due to the occurrence of anearthquake is larger than the immediate increase in risk from an electric signal event.
Another way of examining the relative contribution of the self-exciting and the external excitationRTRT
terms is to compare the contributions of the three terms µ, (1/T ) 0 λS (t)dt and (1/T ) 0 λE (t)dt to theRToverall mean rate (1/T ) 0 λ(t) = NE /T . The results are shown in Table 6. There is substantial variationin the estimates of the background rates and percentages between the different stations. Some of this isprobability due to model instability, but the entries broadly reflect features we have already identified,such as the extensive clustering in the Huailai records.
Overall, both the self-exciting and the electric signal terms make significant contributions to the torisk, with the signals term contributing typically about one-third of the total risk.

6

Supplementary Analyses

The present section describes three additional analyses, all of which were undertaken to investigate aspectsof the stability or interpretation of the previous results. The most important of these, described first,
investigates the possibility that the electric signals are not pre-seismic but post-seismic in character.
This is achieved by reversing the roles of the electric signals and the earthquakes. Then we examine therobustness of the signals effect by allowing a slight over-fitting of the clustering effect in the earthquakes,
and only then testing the effect of adding in the signals as an external term. Finally, we break the period
after the time of the Datong sequences, use the initial period to fit the model, and the later period totest it.

6.1

Analysis of electric signals using earthquakes as explanatory variables

In this analysis we take the same modelling framework as the preceding section, but reverse the roles ofearthquakes and the earthquakes. That is to say, we consider Poisson and self-exciting models for electric
signal events, and then examine the effect of adding in the earthquakes as a possible external exciting
term. The results are summarized in Table 7, which has the same structure and interpretation as Table
4, so that we restrict our remarks to a few key points of interest.
The comparison between the Poisson and self-exciting models shows that the electric signal eventsare even more highly clustered than the earthquakes. The large reductions in AIC values shown in thesecond row block of the table correspond to information gains per event varying from 0.71 for Huailai to
1.37 for Changli.
By contrast, addition of the earthquakes as externally exciting events causes almost no improvementof fit, whether or not the self-exciting terms are included. The third block in Table 7 shows that, evendespite the clustering in both sets of events, for three of the stations AIC selected the Poisson modelin preference to a model with external excitation from the earthquakes. In the final block, when bothself-exciting and external terms are added, the earthquake terms add a significant effect only for Huailai.
8An explanation of the anomalous behaviour of Huailai station is suggested when we plot the conditional intensity function for the electric signals, as in Figure 8(a). From this it is clear that there isan increasing trend in the frequency of electric signals at Huailai. At the same time, the response function for the external exciting term, as shown in Figure 8(c), has such a slow decay that the cumulative
sum builds up over time, thus emulating the trend. In effect, therefore, the externally exciting term for
Hualilai models not the direct predictive effect of an earthquake, but the increasing trend in the numbersof signals detected at the Huailai station. A similar but smaller trend exists also for Langfang station. A
more typical response curve for the self-exciting term for the electric signals is shown in Figure 8(e), for
Qinxian station, where the self-exciting model was selected by AIC.
In no case do we find convincing evidence that the earthquakes trigger the electric signals.

6.2

Testing with an Extended Model for the Self-Exciting (Clustering) Term

As one of the referees pointed out to us, in the present context one of the drawbacks of the Lin-Lin modelis the rapid exponential decay of the self-exciting term, which, in association with AIC selection criterion,
may limit the ability of the self-exciting term to explain subsequent events which are then picked up by
the electric signals. This effect may exaggerate the effectiveness of the electric signals.
As a first step to investigate this effect, we carried out some subsidiary analyses using the ETAS
model to give a better picture of the self-exciting term. The conditional intensity for the ETAS model
has a similar form to (3) except that the decay function is modelled by a power-law term based on the
Omori formula, rather than on the Laguerre polynomial formulation in (4) (see eg Ogata, 1989, 1999) for
further details).
Table 8 shows the results of these analyses. In this table, the results from the ETAS model replacethose from the self-exciting term of the Lin-Lin model, while those from the combined model are left
unchanged from Table 4. Thus the fit from the combined Lin-Lin model is compared to the fit from the
ETAS model. Only for Huailai station does the change affect the character of the results. For all the otherstations, the electric signals still produce a significant improvement in fit.
For Huailai station, it is evident that the ETAS model fits the two large aftershock sequences verymuch better than the self-exciting term from the Lin-Lin model. When compared to the fit from the
combined model, this improvement (which covers half of the total earthquakes in the Huailai catalogue)
swamps any effect from the electric signals. Except for the relatively short period covered by these
two sequences, however, the electric signals appear to precede earthquake events much as for the otherstations. It therefore seems likely that a more complete analysis, which used a combination of the ETAS
and electric signals terms in place of the combined Linlin model, would show that the electric signals
retained some explanatory power.

6.3

Testing with Separate Learning and Evaluation Periods

In Table 9, we summarize the results of fitting the Lin-Lin and Etas models to data from Sanhe station(the one with the longest record) in two parts. The first (training) period, 1982 - end 1989 was used to fit
the Poisson, self-exciting, combined, and ETAS models; the second (evaluation) period (1990 - end 1998)
was used to assess the performance of the model on an independent data set. The combined model gavethe lowest AIC during the training period, and was clearly selected as the best of the four models. The
combined model also gave the highest log-likelihood value in the evaluation period. More than 1/3 of theearthquakes in the Sanhe catalogue (26) occurred during the evaluation period, making the information
gain/earthquake (comparing combined to self-exciting, as in Table 5) about 7.5/26 ≈ 0.28, which is quite
comparable to the overall gains reported in that table. The results for the Huailai station again show thedominance of the clustering term as registered by the ETAS model. When compared to the self-exciting
model, there is still a modest information gain/event.

7

Discussion and conclusions

In this paper, we have used a version of the Hawkes mutually-exciting model to examine the statisticalcorrelations between the electric signal events, which were recorded by 5 stations in the Beijing region,and the intermediate size earthquakes which occurred within a 200-300km radius of each of the stations.
For each station separately, the results indicate that the signals provide weak but nonetheless nontrivial precursory information about the occurrence of magnitude 4 and larger earthquakes in the circular9region around the station. In the converse direction the results are equally unequivocal; in no case do theearthquakes appear to carry significant precursory information about the occurrence of the signal events
recorded at the station.
In response to concerns expressed by the referees and others who have commented on earlier draftsof this paper, we would like to make the following comments.
1. While it may be that the electric signals data analyzed in this study have been collected by relatively
simple equipment, and that the initial data handling requires a degree of subjective intervention
and interpretation, the data analyzed are part of a systematic, standardized process of measurement
and reporting that has been carried out in the Beijing region for the better part of twenty years.
From the point where they have been collected by CAP, the data used are quantitative, objective,
and available for analysis. Within CAP, they form a regular and integral part of the surveillance
and forecasting programme, and have done for many years.
2. During the period of the study, all five stations report very comparable effects, despite local variations in seismicity and geological and man-made environment. Even if the observation period isbroken up into sections, similar effects are observed in each section, as outlined in section 6.3. The
greatest deviations from the common pattern occur for Huailai station. Even here it seems likelythat the discrepancies are caused by the two large aftershock sequences which enter the earthquake
subcatalogue for this station, and that the discrepancies would be much reduced if the aftershock
periods were either better modelled in the combined model, or removed from the analysis.
3. The statistical analysis applied to this data is also crude. Information is lost in converting the signals
to 0-1 data, and only the simple Lin-Lin model has been used in the main analysis. Nevertheless it
is sufficient to show that the precursory effect of the electrical signals persists even after allowance
has been made for earthquake clustering, and that the electric signals are pre- and not post-seismic
in character. Interpretation of the AIC values and associated significance levels is subject to someuncertainty, but in most cases there is a large margin of error. The fact that similar significancelevels are recorded for all stations substantially increases the overall significance.
4. It is clear that in some cases the signals have been contaminated by the occurrence of electric storms,
man-made electrical appliances, etc. Except in the most obvious cases, all such disturbances havebeen incorporated into the assessment of the daily strengths. In such a situation, unless the interfering signals carried predictive information similar to those of the anomalous signals themselves,
one would expect them to act as noise and to decrease the apparent correlations between the signals
and earthquakes. Additional studies by Guan and coworkers suggest that in fact the majority ofthe observed signals are not related to interference from meteorological phenomena or magnetic
storms. For example, the electric signals do not appear to reflect annual fluctuations in rainfall or
temperature (Guan et al, 1996), nor do the large-scale magnetic storms appear to be associated withmajor variations in either the electric signals or the earthquakes (ibid). As for man-made signals,
the two sites most likely to be affected by these, namely Langfang and Huailai, both show increased
frequency of electric signals with time, and are accompanied by rather low values of the predictive
power per signal. These features suggest that, for these stations at least, the earthquake-associatedsignals are indeed diluted by noise signals which build up as urbanization of the area increases. At
least for the observation periods in the present study, however, these effects are not so large as todestroy the correlations with the earthquakes.
5. At present, we have no answer to what may be the most serious problem associated with the present
topic, namely the absence of a plausible physical process that can explain the source and precursory
character of the electric signals. Our view here is that the existence of such a mechanism cannot yet
be ruled out, and that in such a situation the collection and presentation of empirical data, such as
outlined in the present study, is a necessary stage in the scientific process.
6. Although we are aware of the controversy which has surrounded earlier reports of links between
ground electric signals and earthquake occurrence, and the many serious criticisms which have been
directed towards some of these, we believe that the present data and analysis should be given theopportunity to be considered in their own right. At the least the data analyzed here summarize theresults of two decades of preparatory work which is unique of its kind in scale and consistency of
application.

10

In conclusion, we should emphasize that this is only a preliminary analysis of the data on electric
signals held by the China Seismological Bureau Centre of Analysis and Prediction. Its aim was to establishthat there is at least a prima facie case for taking seriously the claims that there may be a connection
between electric signals and earthquakes. At this point, we cannot fully exclude the possibility thatthe positive results obtained in this analysis have their roots in aspects of the data handling, selection
procedures, and analysis, or in some coincidental relation between earthquake occurrence and man-madeor physical noise. Our belief, however, is that this is unlikely to be the case.
We should also emphasize that the present paper does not pretend to set up a quantitative forecasting
scheme. It is based on a retrospective study of existing data, and leaves unresolved many issues relating
to the combination of information from the different stations, and the estimation of locations and sizesof possible future earthquakes, all of which would have to be resolved before a quantitative forecasting
scheme could be put into place. Some of these aspects are pursued further in the companion study by
Zhuang et al (2002), which is currently available as a technical report; some part of its findings are
summarized in Ogata and Zhuang (2001).
At the present time, the electric signals stations around Beijing are in the process of being converted
to digital stations. When this process is completed, it should allow an unequivocal answer to be given toconcerns about reliability and objectivity of the data which have been raised. Should the results confirm
the findings in the present report, the existing data will be a unique source of information about theelectric signals, and their protection becomes a matter of both national and international importance.

Acknowledgements
We would like to acknowledge the critical but nonetheless helpful advice that we have received from theReferees (Yan Kagan and Francesco Mulargia), from Robert Davies, and from other colleagues, all of whichhave contributed substantially to improvements to the original draft of this paper. The work describedin the present study was initialized with assistance of No. 10371012 of the National Natural Science
Foundation of China, then was surported primarily with partial assistance of a subcontract on earthquake
forecasting with the Institute of Geological and Nuclear Sciences on FRST Contract Natural Hazards
of New Zealand and Gran-in-Aid 11680334 for the Scientific Research from the Ministry of Education,
Science, Sport and Culture, Japan. The observation part on the ground ULF electric signals by the fourthauthor has been supported by No. 2001BA601B01-05-03 of the 10th Five-Year-Term Scientific ResearchFunds of China Seismological Bureau. All of these financial assistance are gratefully acknowledged here.
The writing of the paper was assisted by grants from the Centre for Mathematics and its Applicationsin Minneapolis.

References
Akaike, H. (1974), A New Look at the Statistical Model Identification, IEEE Transaction on AutomaticControl, AC-19, 716–723.
Chen, Z., Du, X., Song, R. and Guan H., Electromagnetic Radiation and Earthquakes (in Chinese),(Seismological Press, Beijing, 1998).
Chen, Z. and Xu, C., Features of Electric-Magnetic Radiation Between Some Shallow Earthquakes inNorth China, In Studies on New Method of Short-term Earthquake Prediction, (eds. Department of
Monitoring, State Seismology Bureau), (Seismological Press, Beijing, 1994).
Daley, D. and Vere-Jones, D., An Introduction to the Theory of Point Processes, Vol. 1. (2nd edition).
(Springer, New York, 2002).
Guan, H. and Liu, G. (1995), Electromagnetic Radiation Anomalies Before Moderate and Strong Earthquakes, Acta Seismologica Sinica, 9, 289–299.
Guan, H., Chen, Z. and Yu, S. (1996). Analysis on the Relationship Between Earthquakes and ElectricMagnetic Radiation, Earthquake Research in China, 16, 168–176.
Guan, H., Zhang, H., Liu, Y., Yu, S. and Chen, Z. (1999), Studies on the Relationship Between the
Anomalies of the Ultra-Low Frequency Electric Field and Earthquakes (in Chinese with English
abstract), Earthquake, 19, 142–148.
11

Guan, H., Chen, Z. and Yu, S. (2000), Study on Seismic Precursor Features of Electromagnetic Radiationin the Beijing Region and its Surrounding Area (In Chinese with English abstract), Earthquake,20, 65–70.
Harte, D. and Vere-Jones, D. (2004), The entropy score and its uses in earthquake forecasting, This
volume.
Hawkes, A. G. (1971), Spectra of Some Self-Exciting and Mutually Exciting Point Processes, Biometrika,
58, 83–90.
Kagan, Y. Y. and Knopoff, L. (1977), Earthquake Risk Prediction as a Stochastic Process, Physics ofthe Earth and Planetary Interior, 14, 97–108.
Ma, L. and Vere-Jones, D. (1997), Application of M8 and Lin-Lin Algorithms to New Zealand Earthquake
Data, New Zealand Journal of Geology and Geophysics, 35, 77–89.
Michael, A. (1997), Test Prediction Methods: Earthquake Clustering Versus the Poisson Model, Geophysical Research Letters, 24, 1891–1894.
Ogata, Y., Akaike, H. and Katsura, K. (1982), The Application of Linear Intensity Models to the Investigation of Causal Relations Between a Point Process and Another Stochastic Process, Annals of
the Institute of Statistical Mathematics, 34, 373–387.
Ogata, Y. (1983), Likelihood Analysis of Point Processes and its Application to Seismological data,
Bulletin of the International Statistical Institute, 50, 943–961.
Ogata, Y. (1989), Statistical Model for Standard Seismicity and Detection of Anomalies by Residual
Analysis for Point Process, Tectonophysics, 169, 1–16.
Ogata, Y., (1999), Seismicity Analysis Through Point-Process Modelling: A Review, Pure and AppliedGeophysics, 155, 471–507.
Ogata, Y. and Zhuang, J. (2001), Statistical examination of anomalies for the precursor to earthquakes,
and the multi-element prediction formula: Hazard rate changes of strong earthquakes M ≥ 4.0
around Beijing area based on the ultra-low frequency electric observation (1982-1997), Report of
the Coordinating Committee for Earthquake Prediction, Geographical Survey Institute, Tsukuba,
Japan, 66, 562–570.
Utsu, T. and Ogata, Y. (1997), Statistical analysis of seismicity. In Algorithm for Earthquake Statistics
and Prediction, IAPSPEI Software Library (International Association of Seismologigy and Physics
of the Earth’s Interior. 6, 13—94).
Vere-Jones, D. and Ozaki T. (1982), Some Examples of Statistical Estimation Applied to Earthquake
Data , Annals of the Institute of Statistical Mathematics, 34, 189–207.
Vere-Jones, D. (1998), Probabilities and Information Gain for Earthquake Forecasting, Computational
Seismology, 30, 248–263.
Zhuang, J., Ogata, Y., Vere-Jones, D., Ma, L. and Guan, H. (2002), Statistical confirmation of a relationship between excitation of the utra-low frequency electric signals and magnitude M ≥ 4 earthquakes
in a 300 km radius region around Beijing, Research memorandum 847, Institute of Statistical Mathematics, Tokyo, Japan.

12

Table 1: Observation periods and observed numberS of events at each station
StationHuailaiLangfang
Sanhe
ChangliQingxian

Location
Observation period NEQ NEQ /T NES NES /NEQ40◦ 260 , 115◦ 310 ’87-05-1 ∼ ’98-01-31
59
5.32
239
4.05
39◦ 320 , 116◦ 410 ’82-01-1 ∼ ’96-12-31
77
5.13
366
5.73
39◦ 590 , 117◦ 040 ’82-01-1 ∼ ’98-01-31
75
5.00
266
3.55
39◦ 480 , 119◦ 190 ’90-01-1 ∼ ’96-12-31
22
3.14
42
6.00
38◦ 350 , 116◦ 500 ’82-01-1 ∼ ’96-12-31
60
2.20
271
8.21
NEQ : number of earthquakes; NES : number of electric signal events.
Rates are in events/year.

Table 2: A segment of the records of daily strength from Sanhe station
Date······
1982-04-23
1982-04-24
1982-04-25
1982-04-26
1982-04-27
1982-04-28
1982-04-29
1982-04-30
1982-05-01
1982-05-02
1982-05-03
1982-05-04
1982-05-05
1982-05-06
1982-05-07
1982-05-08
1982-05-09
1982-05-10
1982-05-11
1982-05-12
1982-05-13
1982-05-14······AN S

AEW

AN S + AEW

0
0
0
0
0
434
217
389
401
216
50
389
322
611
16
37
0
12
0
2
1
0

0
0
0
0
0
134
19
167
211
17
0
122
17
235
1
2
0
1
0
0
0
0

0
0
0
0
0
568
236
556
612
233
50
511
339
846
17
39
0
13
0
2
1
0

For the purposes of the present analysis, electric signals are noted as having been present on a given day
if and only if the total strength AN S + AEW exceeds 200. In other cases the signals are noted as having
been absent.

13

Table 3: List of the earthquakes in this analysis. The coding in
the last columkn indicate the stations for which the given event
lay within one or other of the observation circles: H=Huailai, L= Langfang, S = Sanhe, C = Changli, Q = Qinxian. Lower suffix(eg H1 ) refers to smaller event within 200km radius. Upper suffix(1 ) refers to larger event in outer annulus (M > 5 between 200
and 300km from station). Double suffix (11 ) refers to a larger event
within the smaller circle.
Date
1982-01-20
1982-01-26
1982-01-27
1982-02-09
1982-03-08
1982-03-30
1982-07-09
1982-08-09
1982-10-19
1982-12-10
1983-01-31
1983-02-06
1983-04-03
1983-07-21
1983-08-03
1983-08-08
1983-08-21
1983-08-26
1983-09-09
1983-09-15
1983-10-05
1984-01-07
1984-02-27
1984-03-11
1984-03-16
1984-05-17
1984-11-27
1984-11-27
1984-12-13
1985-01-26
1985-04-22
1985-04-22
1985-04-22
1985-05-22
1985-08-11
1985-08-11
1985-08-11
1985-08-23
1985-08-25
1985-10-05
1985-10-17
1985-11-21
1985-11-23
1985-11-30
1986-01-11
1986-02-05
1986-02-15

Time
14:54:13
17:47:57
12:30:42
17:33:37
03:41:59
12:06:33
17:18:24
04:07:11
20:46:01
02:16:46
01:50:03
08:05:47
10:16:30
05:55:55
14:21:15
09:04:58
05:00:31
17:22:18
22:12:45
12:18:35
00:25:17
19:18:22
05:00:08
11:24:07
11:43:54
11:08:26
09:54:20
09:58:20
11:20:25
20:55:53
11:31:33
19:02:37
19:56:59
18:51:44
01:17:16
01:30:32
02:29:15
23:19:31
17:05:19
12:01:58
00:08:25
19:42:22
08:47:14
22:38:00
17:54:35
17:18:36
07:08:47

Latitude37◦ 02037◦ 240
39◦ 45039◦ 38039◦ 52039◦ 48038◦ 070
38◦ 58039◦ 53040◦ 28037◦ 26037◦ 31040◦ 45040◦ 49037◦ 38040◦ 400
39◦ 43039◦ 48040◦ 28037◦ 31039◦ 51039◦ 43039◦ 260
39◦ 38038◦ 28038◦ 58040◦ 30040◦ 29039◦ 39037◦ 22039◦ 45039◦ 440
39◦ 43039◦ 500
39◦ 440
39◦ 45039◦ 460
39◦ 45039◦ 440
39◦ 47039◦ 17040◦ 05038◦ 47037◦ 14037◦ 260
39◦ 45037◦ 450

Longitude117◦ 320114◦ 520
118◦ 460
118◦ 090
118◦ 400
118◦ 300
119◦ 230120◦ 350
118◦ 590
116◦ 330115◦ 010
114◦ 570
114◦ 470114◦ 420
119◦ 060
115◦ 230118◦ 190
117◦ 130116◦ 340
114◦ 160
118◦ 540
118◦ 450118◦ 010
118◦ 220119◦ 050
119◦ 080
114◦ 080
114◦ 160
118◦ 130114◦ 460
118◦ 460
118◦ 470118◦ 460
118◦ 320118◦ 450118◦ 470118◦ 460
118◦ 450118◦ 300
118◦ 270114◦ 390115◦ 500
121◦ 110114◦ 490
114◦ 540
118◦ 260115◦ 140

14

Magnitude
4.2
4.8
4.4
4.7
4.9
4.2
4.1
4.0
5.3
4.9
4.1
4.5
5.1
4.4
4.4
4.1
4.1
4.3
4.0
4.1
4.7
5.2
4.3
4.0
4.0
4.2
4.5
4.1
4.1
4.2
5.0
4.3
4.7
4.7
4.0
4.5
4.0
4.0
4.4
5.0
4.1
4.7
4.5
5.6
4.1
4.0
4.8

RemarkQ1
Q1
S 1 C1L 1 S1 C 1
S 1 C1S 1 C1
C1
C1
L1 S1 C11
H 1 L 1 S1
L 1 Q1L 1 Q1
H11 L1
H1
C1
H 1 L1S 1 C1H 1 L 1 S1 C 1
H 1 L 1 S1
Q1
S 1 C1L1 S1 C11
L 1 S1 C 1
S 1 C1S 1 C1S 1 C1
H1
H1
S 1 C1Q1
L1 S1 C11
S 1 C1S 1 C1S 1 C1S 1 C1S 1 C1S 1 C1S 1 C1S 1 C1H 1 L1 S11 C11 Q1
H 1 L1 Q 1
H 1 L 1 S1
C1
L1 S 1 Q11
L 1 Q1S 1 C1L 1 Q1

Table 3: (continued)Date1986-02-18
1986-07-02
1986-07-02
1986-11-10
1986-12-16
1987-03-21
1987-05-22
1987-06-07
1987-07-16
1987-08-08
1987-08-08
1987-11-11
1988-05-09
1988-07-23
1988-07-25
1988-07-26
1988-07-30
1988-08-03
1988-10-26
1989-06-18
1989-07-23
1989-10-05
1989-10-18
1989-10-18
1989-10-19
1989-10-19
1989-10-19
1989-10-19
1989-10-19
1989-10-19
1989-10-19
1989-10-19
1989-10-20
1989-10-20
1989-10-23
1989-10-29
1989-12-09
1989-12-13
1989-12-14
1989-12-25
1989-12-31
1990-03-16
1990-05-23
1990-07-21
1990-07-23
1990-08-03
1990-09-22
1990-12-24
1991-01-29
1991-03-26
1991-03-26
1991-04-01
1991-05-07
1991-05-29

Time
06:07:03
11:56:58
11:58:29
16:58:05
07:25:56
07:00:23
11:50:42
12:26:23
02:16:28
07:35:06
17:35:08
21:18:39
06:22:01
13:51:43
01:47:39
22:46:39
18:44:32
17:44:11
01:29:44
14:50:59
03:15:09
11:55:35
22:57:23
23:29:59
00:52:55
01:01:34
01:09:43
01:11:17
01:26:26
02:20:46
05:02:03
18:29:03
01:56:48
19:41:41
21:19:32
10:22:43
07:04:51
09:10:11
12:13:01
04:26:23
16:24:48
15:08:39
14:13:01
08:41:51
16:41:32
18:05:44
11:02:19
13:31:03
06:28:04
02:02:38
02:07:27
05:35:23
00:25:01
19:02:10

Latitude37◦ 48038◦ 23038◦ 24040◦ 03037◦ 160
38◦ 17040◦ 00039◦ 32039◦ 47039◦ 200
39◦ 20040◦ 17039◦ 18040◦ 05039◦ 31039◦ 340
39◦ 460
39◦ 360
39◦ 52039◦ 40037◦ 22039◦ 370
39◦ 570
39◦ 570
39◦ 570
39◦ 570
39◦ 540
39◦ 58040◦ 000
39◦ 59039◦ 55039◦ 55039◦ 560
39◦ 59039◦ 55039◦ 560
39◦ 53039◦ 56037◦ 36040◦ 200
39◦ 58039◦ 46040◦ 13040◦ 35039◦ 45037◦ 53040◦ 05037◦ 54038◦ 28039◦ 58039◦ 59039◦ 540
39◦ 48039◦ 430

Longitude
115◦ 110120◦ 290120◦ 300
116◦ 430114◦ 450114◦ 150
116◦ 470118◦ 050
118◦ 380117◦ 530117◦ 540114◦ 480118◦ 060
114◦ 130118◦ 060
118◦ 060
118◦ 430
118◦ 390118◦ 290118◦ 190
115◦ 020
118◦ 150
113◦ 500113◦ 500113◦ 480113◦ 490
113◦ 460
113◦ 490
113◦ 510
113◦ 520
113◦ 490
113◦ 440113◦ 480113◦ 540113◦ 490
113◦ 490
113◦ 490
113◦ 530115◦ 170118◦ 570
113◦ 510
118◦ 270116◦ 280115◦ 500
118◦ 290115◦ 020
116◦ 220115◦ 010
112◦ 320113◦ 510
113◦ 500113◦ 490
118◦ 420
118◦ 180

15

Magnitude
4.5
4.6
4.3
4.7
4.2
4.5
4.0
4.5
4.5
4.7
4.1
4.7
4.2
5.0
4.4
4.9
4.3
4.6
4.3
4.3
4.3
4.2
5.7
4.2
4.2
6.1
4.0
5.1
4.4
5.6
4.4
5.1
4.1
4.1
5.2
4.1
4.2
4.1
4.8
4.7
4.0
4.2
4.3
5.0
4.9
4.2
4.5
4.1
5.5
6.1
4.3
4.0
4.3
5.2

Remark
L 1 Q1C1
C1
H 1 L 1 S1Q1
L 1 Q1H 1 L 1 S1
L 1 S1 C 1
S 1 C1L 1 S1 C 1
L 1 S1 C 1
H 1 L1L 1 S1 C 1
H11 L11 S 1 Q1
L 1 S1 C 1
L 1 S1 C 1
S 1 C1S 1 C1S 1 C1S 1 C1Q1
S 1 C1
H11 L1 S 1 Q1
H1
H1
H11 L1 S 1 Q1
H1
H11 L1 S 1 Q1
H1
H11 L1 S 1 Q1
H1
H11 L1 S 1 Q1
H1
H1
H11 L1 S 1 Q1
H1
H1
H1
L 1 Q1C1
H1
S 1 C1H 1 L 1 S1
H11 L11 S 1 Q1
S 1 C1L 1 Q1H 1 L 1 S1
L 1 Q1H1
1 1 1 1H1 L S Q
H1H1
S 1 C1H 1 L1 S11 C11 Q1

Table 3: (continued)Date
1991-05-30
1991-07-11
1991-07-27
1991-08-21
1991-08-22
1991-09-02
1991-09-20
1991-09-28
1991-09-28
1991-10-05
1991-10-07
1991-10-17
1992-02-14
1992-07-22
1993-06-14
1993-08-30
1993-08-30
1993-09-27
1993-11-18
1994-10-04
1994-12-23
1995-06-27
1995-07-20
1995-10-06
1995-11-13
1996-04-08
1996-12-16
1996-12-16
1997-04-12
1997-05-25
1998-01-10
1998-01-10
1998-01-10
1998-01-10
1998-01-10
1998-01-10
1998-01-10
1998-01-10
1998-01-11
1998-01-11
1998-01-12
1998-01-12
1998-01-14
1998-01-14
1998-01-17
1998-01-18
1998-01-18
1998-01-22
1998-01-27
1998-02-13
1998-04-14
1998-04-14
1998-04-16
1998-06-02

Time
07:06:55
19:05:05
17:54:47
05:28:28
06:23:37
04:17:42
17:52:06
08:37:47
18:48:06
05:44:39
07:25:21
02:19:42
17:08:43
05:43:02
03:01:00
07:27:59
08:13:38
15:44:16
07:05:09
15:54:17
13:13:40
09:39:51
20:51:23
06:26:53
14:33:19
00:39:28
05:36:33
09:52:24
15:05:02
14:59:08
11:50:39
11:59:17
12:09:58
13:03:59
15:38:18
20:01:55
21:50:32
23:33:20
02:37:47
11:31:43
02:41:54
18:32:29
01:17:26
11:09:57
10:41:18
04:07:25
06:17:54
12:11:51
07:29:07
07:05:05
10:47:48
10:48:06
07:13:44
09:32:04

Latitude39◦ 41039◦ 41039◦ 51037◦ 20037◦ 28038◦ 49039◦ 20040◦ 05040◦ 04039◦ 18037◦ 560
39◦ 440
39◦ 460
39◦ 17037◦ 25039◦ 540
39◦ 540
39◦ 39039◦ 360
39◦ 44040◦ 28038◦ 14040◦ 160
39◦ 400
39◦ 22039◦ 51040◦ 10040◦ 060
38◦ 17040◦ 42041◦ 06041◦ 20041◦ 06041◦ 05041◦ 09041◦ 07041◦ 11041◦ 04041◦ 05041◦ 05041◦ 08041◦ 12041◦ 04041◦ 11041◦ 06041◦ 09041◦ 07041◦ 08041◦ 02041◦ 03039◦ 41039◦ 41041◦ 05041◦ 110Longitude118◦ 160
118◦ 230118◦ 480114◦ 420114◦ 440119◦ 590
114◦ 050
117◦ 030
117◦ 040
117◦ 570
115◦ 040
118◦ 250
118◦ 260117◦ 560
119◦ 560
113◦ 490
113◦ 490
118◦ 420
117◦ 270
118◦ 260115◦ 330119◦ 280115◦ 260118◦ 200113◦ 120
118◦ 440116◦ 300
116◦ 350120◦ 290114◦ 520114◦ 180114◦ 320114◦ 280114◦ 300
114◦ 320114◦ 260114◦ 260114◦ 250114◦ 250114◦ 260114◦ 270114◦ 280114◦ 270114◦ 310114◦ 260114◦ 220114◦ 280114◦ 260114◦ 290114◦ 270
118◦ 280
118◦ 280114◦ 270114◦ 240

16

Magnitude
5.6
4.3
4.4
4.3
4.1
4.1
4.1
4.0
4.0
4.0
4.1
4.3
4.4
4.9
4.1
4.1
4.0
4.0
4.4
4.0
4.3
4.0
4.1
5.4
4.5
4.0
4.5
4.0
4.3
4.7
6.2
4.0
4.5
4.6
4.2
4.0
4.3
4.0
4.6
4.4
4.3
4.3
4.2
4.3
4.3
4.6
4.3
4.4
4.2
4.0
5.0
4.4
4.0
4.8

Remark
H 1 L1 S11 C11 Q1
S 1 C1S 1 C1Q1
Q1
C1
H 1 L1H 1 L 1 S1
H 1 L 1 S1
L 1 S1 C 1
L 1 Q1S 1 C1S 1 C1L 1 S1 C 1
C1H1
H1
S 1 C1L 1 S1 C 1
S 1 C1H 1 L1C1
H 1 L 1 S1
H 1 L1 S11 C11 Q1
H1
S 1 C1H 1 L 1 S1
H 1 L 1 S1
C1
H1
H11 L1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H1
H 1 L1 S11 C11 Q1
S 1 C1
H1
H1

Table 4: Results from fitting models to the earthquake data
StationNP

HuailaiLangfang
Sanhe
ChangliQingxian
HuailaiLangfang
Sanhe
ChangliQingxian
HuailaiLangfang
Sanhe
ChangliQingxian
HuailaiLangfang
Sanhe
ChangliQingxian

NS NE k plog L
Trend model2
2 -33.66
1
1 -50.77
1
1 -56.67
1
1 -25.29
1
1 -54.52
Self-exciting model
1
1
3
89.10
1
1
3 -22.77
1
1
3 -39.98
1
1
3 -23.11
1
1
3 -42.23
Externally excited model
1
2
4
15.68
1
2
4 -30.14
1
2
4 -32.20
1
1
3 -20.13
1
2
4 -26.08
Combined model
1
1
2
6 100.86
1
1
2
6
-8.70
1
1
2
6 -15.83
1
1
1
5 -17.67
1
1
2
6 -22.65

AIC
71.31
103.54115.3452.58
111.06-170.22
51.55
83.96
52.22
90.47
-23.35
68.28
72.46
46.26
60.17
-189.72
29.40
43.66
45.34
57.31

∆kp

log L/L0 ∆AIC
vs. trend model
1
122.76 241.53
2
28.00 -51.99
2
16.69 -31.38
2
2.18
-0.36
2
12.28 -20.59
vs. tend model
2
49.34 -94.66
3
20.63 -35.26
334.47 -42.88
2
5.16
-6.32
3
28.44 -50.99
vs. Self-exciting model
3
11.76 -19.50
3
14.07 -22.15
3
24.15 -40.30
2
5.44
-6.88
3
19.58 -33.16

NS : order of the Laguerre polynomial used for the self-exciting term; NE : order of the Laguerre polynomial used forthe external excitation term; kp : total number of fitted parameters in the model; log L: maximized log likelihood;
AIC: as defined by (8); ∆p: difference between the numbers of parameters in the two models for comparison;log L/L0 : log likelihood ratio of the two models for comparison; ∆AIC: difference in AIC values between the twomodels for comparison;

Table 3: (continued)Date1998-06-02
1998-07-14
1998-07-27
1998-07-27
1998-07-27
1998-08-13
1998-08-15

Time
10:47:29
18:16:08
09:05:16
09:16:56
11:17:48
18:21:52
18:11:22

Latitude41◦ 03041◦ 06041◦ 12041◦ 06041◦ 10041◦ 01038◦ 380

Longitude114◦ 300
114◦ 270114◦ 290114◦ 270114◦ 290114◦ 350
119◦ 140

17

Magnitude
4.5
4.4
4.8
5.0
4.7
4.7
4.4

RemarkH1
H1
H1
H11 L1
H1
H1
S 1 C1

Table 5: Information gains per earthquake and per signal for each station
StationInfogain/earthquake Infogain/earthquake
Infogain/signal(cluster v Poisson)(combined v cluster) (combined v cluster)
Huailai
2.08
0.20
0.049
Langfang
0.36
0.18
0.038
Sanhe
0.22
0.28
0.091
Changli
0.10
0.25
0.130
Qingxian
0.20
0.33
0.072
The information gain means here the log-likelihood ratio, normalized by
dividing by the number of earthquake events or the number of signal events.

Table 6: Comparison of intensities in the combined model.
RT
Stationµ
Background Clustering Signals
λE dt/NS
0
Huailai
0.3104
20.6%
50.8%
28.6%
0.071
Langfang 0.7601
57.2%
14.9%
27.9%
0.056
Sanhe
0.6814
53.4%
10.5%
36.2%
0.102
Changli
0.4143
48.1%
29.7%
22.2%
0.116
Qingxian 0.4084
37.3%
15.7%
47.1%
0.104
µ is the estimated background rate in events per 100 days; the next three columns give percentages of
the total rate due to each term in the intensity; the final column is a rough indication of the number ofpredicted earthquakes per electric signal.

Table 7: Results from fitting models to the electric signals data
St.
HuailaiLangfang
Sanhe
ChangliQingxian
HuailaiLangfang
Sanhe
ChangliQingxian
HuailaiLangfang
Sanhe
ChangliQingxian
HuailaiLangfang
Sanhe
ChangliQingxianNP

NS NE k plog L
Poisson model2
2 195.23
22 479.90
1
1 135.76
1
1 -21.12
1
1 162.32
Self-exciting model
2
1
4 378.04
2
1
4 741.20
12
4 384.47
1
1
3
28.45
1
2
4 433.99
Externally excited model
1
3
5 221.16
2
1
4 483.47
1
3
5 137.58
1
1
3 -21.12
1
9 11 198.01
Combined model
12
2
7 391.29
2
1
1
6 741.20
12
1
6 384.65
1
1
1
5
30.34
1
2
2
7 434.56

AIC-386.46-955.79-269.5244.25-322.65-748.08-1474.40-760.94-50.90-859.98-432.32-958.95-267.1648.25-334.01-768.59-1470.40-757.30-50.68-857.13∆plog L/L0 ∆AIC
vs. Trend model2
182.81 -361.62
2
261.30 -518.61
3
248.71 -491.42
2
49.57-95.15
3
271.67 -537.33
vs. Trend model
3
25.93-45.86
2
3.57
-3.16
4
1.82
2.36
2
0.00
4.00
10
35.69
-11.36
vs. Self-exciting model
3
13.25-20.51
2
0.00
4.00
2
0.18
3.64
2
1.89
0.22
3
0.57
2.85

NS : order of the Laguerre polynomial used for the clustering term; NE : order of the Laguerre polynomial used
for the external excitation term; kp : total number of fitted parameters in the model; log L: maximized logarithm
likelihood; AIC: as defined by (8); ∆p: difference between the numbers of parameters in the two models; log L/L 0 :logarithm likelihood ratio of the two models; ∆AIC: difference in AIC values between the two models.

18

Table 8: Fitting results from the ETAS model
StationHuailaiLangfang
ShanheChangliQingxianµ
0.33461.008.8576
0.74290.8161

K
0.002340
0.006772
0.020220.018070.005854c0.0001442
0.004906
0.0007426
0.01675
0.02902

α
2.608
1.012
0.34380.00000.2707

p
1.088
1.200
1.094
1.423
2.029

log L
120.2
-16.73-35.64-28.10-42.32AIC-230.4
43.45
81.28
66.20
94.63

AICl represent the AIC-values for the corresponding best combined model in Table 4.

Table 9: Training and evaluation

Model
Combined
Self-exciting
ETAS
PoissonModel
Combined
Self-exciting
ETAS
PoissonSanhe
Training (NEQ = 49)
1982.1.1 – 1989.12.31
log L
AIC
0.87
10.26
-11.07
28.15
-8.99
27.98
-23.6749.34
Huailai
Training (NEQ = 23)
1982.1.1 – 1990.12.31
log L
AIC
43.18
-74.3538.16-70.3256.90-103.80-10.5923.18

19

Prediction (NEQ = 26)
1990.1.1 – 1998.1.31
log L-23.67-31.20-29.08-36.01
Prediction(NEQ = 36)
1991.1.1 – 1998.1.31
log L
48.18
48.54
60.37
-24.95AICl-189.72
29.40
43.66
59.62
57.31

43
42
41
40

S

C

Beijing

39

L

37

38

Q

36

Longitude

H

112

114

116

118

120

122

LatitudeFigure 1: Distribution of earthquakes and electric signals observation stations around Beijing. Observation
stations are denoted by black squares, with station initial inside. The two circles represent 200- and 300km
circles around Lanfang station. Earthquakes are represented by small circles. The diameter corresponds
to size of event, from ML = 4 for smallest to ML = 6.7 for largest. Dark circles denote events whichfell within the observation region for at least one station; light circles denote other events with similar
magnitudes. The Tangshan area is to the East, adjacent to the Changli station. The Datong and Zhangbei
aftershock sequences are the clusters to the S.W and N.W. of Huailai station, respectively.

20

(a)

(b)

(c)

Figure 2: (a) A general circuit diagram for the observation equipment; (b) A photograph of the electrodes
and the DJ-1 recorder. The body of the sensor of the electrodes is made of high quality Cr18Ni9C stainlesssteel. The water-proof pipe (white) will cover the whole screened transmission cables (black) when theelectrodes are installed; (c) The DJ-1 recorder.

21

1
0.5
0.01

Response
0.05
0.2
0.05

0.2
1
5
20
Periods in seconds

Figure 3: The observation system response function.

Figure 4: A photocopy of an original record. The fluctuations in the upper part of the chart are typical of
the observed anomalous signals. The large pulse in the lower half of the chart is an example of a clearly
distinguishable interference pulse, not included in the computation of the daily signal strength.

22

3 4 5 6
3 4 5 6

Magnitude

3 4 5 63 4 5 63 4 5 6(a) Huailai(b) Langfang(c) Sanhe(d) Changli(e) Qingxian

1982 1984 1986 1988 1990 1992 1994 1996 1998

Time (Year)
Figure 5: Electric signal data and earthquake data used in this analysis. Electric signals (ie days where 1
is recorded) are represented by heavy vertical lines reaching to M=4. Earthquakes are represented by thin
vertical lines, the height corresponding to ML on the left hand scale. Cumulative numbers of earthquakes
are represented by the thick grey lines, scaled by the dashed horizontal lines from 0 to 80. Periods withno events recorded correspond to periods for which the stations were not operational.

23

(a)

(b)

Figure 6: A photocopy of a record with noises (the upper panel) from ground strong lightnings and (thebottom panel) from electricity leaked by power lines. In the upper panel, the short vertical bars, which
form up a line intersecting the horizontal line, are the timing marks at every minute, and irregular pulses
caused by strong ground lightnings are marked by arrows. In bottom, the part in black represent theperiod of power leakage.

24

1.2
0.9
0.6
0.3
0

Events/day

(a)

1982

1984

1986

1988

1990

1992

1994

1996

1998

0.3

0.6

0.9

(b)

0

Events/day

1.2

Time (Year)

1989

1990

1991

1.2
0.8
0.0

Events/100days

(d)

0.4

10 20 30 40 50(c)0Events/100days

Time (Year)

0.00

0.05

0.10

0.15

0.00

Time (/100 days)

0.05

0.10

0.15

Time (/100 days)

Figure 7: Fitting results for the earthquake data around the Sanhe station using the ULF electric signal
events as explanatory factors. Top: Conditional intensity function throughout the whole time period;
Middle: Conditional in a short period; Bottom left: impulse response from self excitation; Bottom right:
impulse response from mutual excitation.

25

1.7
0.7 1 1.3
0 0.3

Events/day

(a)

1987

1989

1991

1993

1995

1997

0.12
0.08

(c)

0.00

0.04

Events/100days

5

10

15

(b)

0

Events/100days

Time (Year)

0.00

0.05

0.10

0.15

0

20

30

40

50

0.1

0.3

(d)

0

Events/day

10

Time (/100 days)

0.5Time (/100 days)

1982

1984

1986

1988

1990

1992

1994

1996

5

10

15

(e)

0

Events/100days

Time (Year)

0.00

0.05

0.10

0.15

Time (/100 days)

Figure 8: Fitting results for the ULF electric signal events from the Huailai and Qianxian station using theearthquake data as explanatory factors. (a) Conditional intensity function at Huailai station; (b) impulse
response from self excitation at Huailai station; (c) impulse response from external excitation at Huailai
station; (d) conditional intensity function at Qinxian station; (e) impulse response from self excitation at
Huailai station (no response from external excitation because the model with only self-exciting terms is
selected).

26

REGISTRATION
of

Research Memorandum

Research Memorandum NO.

916 ;(month) (day) (year)
06 , 14 , 2004
Received onbyCenter for Information on Statistical Sciences
the Institute of Statistical Mathematics
Applicant: Zhuang Jiancang

phone: 03-5421-8780, email: zhuangjc@ism.ac.jp·····················································································

Title:

Preliminary Analysis of Observations on the Ultra-Low Frequency Electric
Field in a Region around Beijing
Author(Affiliation):

Zhuang, Jiancang (Institute of Statistical Mathematics, 106-8569, Tokyo,
Japan);
Vere-Jones, David (School of Mathematical and Computing Sciences, Victoria
University of Wellington, PO Box 600, Wellington, New Zealand);
Guan, Huaping (Center for Analysis and Prediction, China Seismological Bureau, 100036, Beijing, China);
Ogata, Yosihiko (Institute of Statistical Mathematics, 106-8569, Tokyo, Japan);
Ma, Li (Center for Analysis and Prediction, China Seismological Bureau,
100036, Beijing, China)
Key words:

ultra-low frequency electric signal; earthquake; precursor; self-exciting andmutually exciting model
Abstract:

This paper presents a preliminary analysis of observations on ultra-low
frequency ground electric signals from stations operated by the China Seismological Bureau over the last 20 years. A brief description of the instrumentation, operating procedures and data processing procedures is given. The data
analysed consists of estimates of the total strengths (cumulated amplitudes) of
the electric signals during 24 hour periods. The thresholds are set low enoughso that on most days a zero observation is returned. Non-zero observations
are related to electric and magnetic storms, occasional man-made electricaleffects, and, apparently, some pre-, co-, or post-seismic signals. The main
— to be continued —Registration page 1

— continued —purpose of the analysis is to investigate the extent that the electric signals
can be considered as pre-seismic in character. For this purpose the electric
signals from each of five stations are jointly analyzed with the catalogue of
local earthquakes within circular regions around the selected stations. A version of Ogata’s Lin-Lin algorithm is used to estimate and test the existence ofa pre-seismic signal. This model allows the effect of the electric signals to betested, even after allowing for the effects of earthquake clustering. It is foundthat, although the largest single effect influencing earthquake occurrence isthe clustering tendency, there remains a significant preseismic component
from the electrical signals. Additional tests show that the apparent effect is
not post-seismic in character, and persists even under variations of the modeland the time periods used in the analysis. Samples of the data are presented,
and the full data sets have been made available on local websites.

Registration page 2