CC BY
9
RUSSIAN JOURNAL OF EARTH SCIENCES, VOL. 12, ES5003, doi:10.2205/2012ES000522, 2012
Method of spectral-time analysis for recognition of anomalies in time series with Raleigh- and tsunami-wave disturbances in signals from hydrostatic pressure sensors of ocean bottom seismic stations
V. G. Getmanov1
Received 26 July 2012; accepted 29 July 2012; published 12 September 2012.
The article deals with the problem of recognition of anomalies in time series with Raleigh- and tsunami-wave disturbances in signals from hydrostatic pressure sensors (HPS) of ocean bottom seismic stations in the framework of tsunami warning problem. The proposed method of spectral-time analysis (STAN) of signals from the sensors was based on computing the evaluations of functions of frequency-time distribution, decision-making procedures and non-linear filtering for the above problem. The developed STAN method was applied to recognize time intervals with Rayleigh and tsunami-wave disturbances in HPS signals. The proposed STAN method is quite universal and can be used to solve problems of recognition of anomalies in time series of geophysical data of different nature. KEYWORDS:
Tsunami; hydrostatic pressure sensors; recognition; anomalies; Raleigh disturbances; spectral-time analysis; functions of frequency-time distribution; decision-making procedures.
Citation: Getmanov, V. G. (2012), Method of spectral-time analysis for recognition of anomalies in time series with Raleigh- and tsunami-wave disturbances in signals from hydrostatic pressure sensors of ocean bottom seismic stations, Russ. J. Earth. Sci., 12, ES5003, doi:10.2205/2012ES000522.
1. Introduction
Recognition of time series in geophysical data, resulting from various geophysical processes, is one of the most important areas in modern physics of the Earth. Typically, activation of geophysical processes manifests itself in the formation of anomalies in time series of geophysical data. Thus, short-time sections in seismograms, corresponding to the processes of earthquakes, are frequently seen as abnormal, and time series of the Earth’s magnetic field magnetograms, correlating to the processes of magnetic storms, also can be interpreted as anomalies.
The general principles of the task of recognition of anomalies in time series of geophysical data refer to [Gvishiani and Dubois, 2002]. In the works [Agayan et al., 2005; Bogout-dinov et al., 2010; Gvishiani et al., 2008] the methods of recognition of anomalies in time series of geophysical data (magnetic in general) were developed. The methods were based on a nonstandard approach, which combines a significant expansion of the initial requirements for the objects of recognition and several provisions of the problem in terms of fuzzy mathematics.
1Geophysical Center RAS, Moscow, Russia
Copyright 2012 by the Geophysical Center RAS. http://elpub.wdcb.ru/journals/rjes/doi/2012ES000522.html
The article confirms the results, obtained in the above works and describes the recognition task from the standpoint of spectral-time analysis (STAN) [Boashash, 2003; Cohen, 1989; Hlawatch and Auger, 2008], which is an effective tool for the study of non-stationary vibration signals. It should be noted that different STAN methods have been successfully used in geophysics, for example, [Dzienovski et al., 1969; Kedrov and Kedrov, 2006; Kedrov et al., 1998; Lander et al., 1973]. The present article is an attempt to combine STAN with decision-making procedures and the subsequent non-linear filtering. The article’s results are aimed at solving the problem of automatic detection of geophysical signals and are close to the content of the work [Baranov, 2007], which describes an alternative technology based on wavelet transforms.
The algorithms of recognition, described in the article, make it possible to compare the results of recognition of anomalies based on fuzzy mathematics, wavelet transforms and the proposed STAN method.
Recognition of time series with Raleigh- and tsunami-wave disturbances (further, LR and TW-disturbances) in signals from hydrostatic pressure sensors (HPS) of ocean bottom seismic stations (BSS), containing signals from earthquakes and tsunamis, is an actual task of geophysics. Raleigh waves are formed in the foci of earthquakes, they are distributed on the bottom surface; HPS signals with such components are early precursors of tsunami. Tsunami-wave disturbances,
Figure 1. Filtered HPS signal for 6, 503, 000 < i < 7, 752, 000.
detected in HPS signals from BSS, which are located far away from the shore, may be taken as direct precursors of a tsunami.
BSS with HPS have many advantages over traditional bottom stations, which are equipped with seismographs [Bashilov et al., 2008]. Certain information on the characteristics of signals and the construction of the BSS with HPS are given in [Dykhan et al., 1981; Kulikov and Gonzales, 1995; Levin and Nosov, 2005]. Detailed information about the construction of American BSS with HPS of DART-2 type (Deep-ocean Assesment and Reporting of Tsunami) is placed on the site [www.ndbc.noaa.gov]. BSS, integrated into the global system [www.ndbc.noaa.gov], can provide an effective tsunami warning and make a significant contribution to the study of seismicity of the Earth. These statements are based on the fact that ~ 80% of all earthquakes occur beneath the bottom of the oceans and seas, and a network of exclusive land-based seismic stations cannot record earthquakes without gaps.
The present article describes a system of sequential algorithms, providing recognition of time series with Rayleigh- and tsunami-wave disturbances in the records of BSS signal. The first STARTS-1 algorithm (Spectral-Time Analysis of Rayleigh-waves and Tsunami-waves Signals) calculates the evaluation of time-frequency distributions (TFD) functions and implements the decision making procedure (d.m.p.) for recognition by comparing the evaluations and reference TFD functions. The second STARTS-2 algorithm produces a non-linear filtering of the d.m.p. results to improve the accuracy of estimating the boundaries of the anomalous areas and to reduce the probabilities of false alarms and omissions.
Some results of this work are described in [Getmanov et al., 2011]. The material of this work in a certain way correlates with [Chebrov and Gusev, 2010; Poplavsky et al., 1988].
2. Signals of Hydrostatic Pressure Sensors, Raleigh- and Tsunami-Wave Disturbances
The records of signals of the National Geophysical Data Center (NGDC-USA) were used, available in the Internet [www.ngdc.noaa.gov]. For the calculations the record of HPS signal from the DART-2 buoy 46,419, located in the northern part of the U.S. West Coast (coordinates 48.4785°N, 129.3593°W) at a depth of ~ 4138 m, was used. The period of record was ~ 3.5 years (2006.10.252010.03.11), quantization step T = 15 s, recording capacity ~ 7, 000, 000 numbers. The measurement accuracy of HPS from DART-2 has adopted a value of about 1 mm of water column [www.paroscientific.com].
Figure 1 shows the implementation of the oscillation signal y(Ti), obtained after filtration of the additive low-frequency tidal components in the HPS signal for the points 6, 503, 000 < i < 7, 752, 000; units of measurement - meters of water column. Observations dates: 2009.08.25,
21 h. 20 min. - 2010.03.30, 17 h. 25 min. The number of points for the signal in Figure 1 was 1,249,000, the period of observations ~ 216.8 days. It was clear that the values of noise components of pressure in the HPS signal were equal, on average, to 3 ^ 5 mm. Several time series of the observed signal y(Ti) corresponded to different variants of disturbances.
Let us consider the variants of disturbances, observed in HPS signals. Let us use the information on the physics of signals [Meining et al., 2005; Tompson, 2009] and online materials [www.ngdc.noaa.gov]. The implemented version of disturbance in the first approximation can relate to a certain time in the signal area.
Figure 2. Filtered HPS signal for the points 6, 703,100 < i < 6, 709, 900.
2.1. Disturbance for variant 1 related to the Raleigh (LR) wave processes. LR-disturbances in the signal occur due to submarine earthquakes (sometimes terrestrial earthquakes) and Rayleigh surface transverse waves that propagate on the bottom surface at a speed of 2000-3000 m s-1. LR-disturbances in the observed signal are caused by pressure fluctuations of the aquatic environment due to seismic effects in the vicinity of BSS.
Figure 2 shows the signal y(Ti)for the points 6, 703,100 < i < 6, 709, 900, which allows to see the LR-disturbance at a larger scale, registered by BSS as oscillation impulses in the area of the points 6, 704,100 < i < 6, 704, 250.
Figure 3 shows the same signal y(Ti) for 6, 704,165 < i < 6, 704, 340 at a larger scale. An oscillatory impulse of the LR-disturbance can be seen, related to the average point
i ~ 6, 704, 200. The LR-disturbance consists of a number of high-frequency oscillations with increasing and then decreasing amplitudes; LR-disturbance, in this case, takes time, approximately equal to 375 s.
2.2. The disturbance for variant 2 is related to tsunami waves-TW-random oscillations, with durations, which
are equal, on average, from 100 ^ 300 to 3000 ^ 5000 points ((1500 ^ 4500) — (4500 ^ 7500) s). In Figure 2 the TW-disturbance corresponds to the section with the points « 6, 706, 750 < i < 6, 708, 000. A TW-disturbance is caused by wave oscillations at the vicinity of BSS. The reason of oscillations relates to seismic vertical oscillatory motions of the ground near the earthquake source, which are then transferred into the aqueous medium and distributed in the water at a rate of 200 m s-1. As a general rule, sites with LR- and TW-disturbances in the observed signals are separated in time due to their different propagation velocities. Figure 4 on a larger scale shows a TW-disturbance, which begins at the points i « 6, 706, 750. We see that a TW-disturbance corresponds to low-frequency oscillations. When considering at a fine time scale, a TW-disturbance is characterized by a sharp leading edge and a flat trailing edge.
2.3. In case of the variant 3, the disturbance is of a double nature. CGM-oscillations-Continuous Ground Motion - are caused by micro-seismic broadband noise from the bottom surface. TWGM (TW + Ground Motion) oscillations are caused by a combination of a tsunami wave
Figure 3. Signal with LR-disturbances at points 6, 704,165 < i < 6, 704, 340.
Figure 4. Signal with TW-disturbance, starts at the point i « 6, 706, 750.
and random noise of seismic effects on the bottom surface. The amplitudes of oscillations for TWGM-disturbances can vary widely; TWGM-disturbances are implemented at sufficiently large time intervals. In Figure 2 the intervals
6, 703,100 < i < 6, 704,100; 6, 704, 250 < i < 6, 706, 750 correspond to CGM-disturbances, TWGM-disturbances correspond to points i > 6, 708, 000. In a number of cases, TWGM-disturbances are a continuation of TW-disturbances. CGM and TWGM-disturbances, in general, are broadband random oscillations.
It must be kept in mind that a real picture, consisting of a set of options for disturbances in the observed signals from the BSS can be quite complicated, for example, when registering vibrations from several simultaneous earthquakes, etc.
LR and TW-oscillations differ in their frequency characteristics from CGM and TWGM-oscillations; they take significantly less time than CGM and TWGM-oscillations. Therefore, in relation to the entire signal, time series with LR, and TW-disturbances, for further convenience of terminology will be interpreted as anomalies.
In this paper we consider the problem of recognition of time series with anomalies of Rayleigh and tsunami-wave disturbances in HPS signals. These signals are non-stationary with time-varying spectral characteristics. For their recognition we use the STAN algorithm combined with the d.m.p. algorithms and nonlinear filtering.
erates time series of geophysical data, is defined on a time interval of observation t0 < t < tf. We also assume that anomalous disturbances in the signal are located in the sections, which are given by initial and final moments of time (t1i, t2i), I = 1, 2,... ,L - the number of anomalies, Figure 5.
The moments of time t, not located in the anomalous sections, will be considered as belonging to non-anomalous sections of the records. Thus we can assume that the initial time interval t0 < t < tf consists of sets of anomalous and non-anomalous sections.
For the signal x(t) we make a provision for its frequency representation based on the Fourier transform C(ju). The basis of the STAN algorithm consists of the so-called functions of time-frequency distributions (TFD) P(u,t), defined in the simplest case, in the rectangle u e Q, t e T, T = {t : (t0 < t < tf)}, Q = : (u>0 < u < Uf)}. The
TFD functions are functions of two variables - time and frequency. The physical meaning of TFD functions follows from the relations P(u,t)dtdu = dE
/ / P (u,t)dtd(j = E
which determine a value dE of a signal energy attributable to a time interval (t,t + dt) and frequency range (u,u + du), and the total energy of signal E . It should be noted that the TFD function and value of the instantaneous signal power calculated in the time and frequency domains are connected by the obvious relations
3. STAN and Recognition of Anomalies in Time Series of HPS Signals
P(u, t)dt = IC(jw)
We present the necessary information from the spectraltime analysis (STAN) and formulate the general problem of recognition. Let us consider the model descriptions of the initial signal and anomalies for the continuous case.
We assume that the initial HPS signal x(t), which gen-
Figure 5. Location of anomalies.
:,v
2
/ P(u, t)du = |*(i)|
J
Let us introduce, using the time points tn, n = 1,... , n0, tno = tf small local time intervals dt„ = {t : (tn — 1 < t < tn)} for t0 < t < tf. Let us determine in local intervals u e Q, t e dt„ the local TFD functions Pn(u,t), which are equal to zero outside the local intervals - Pn(u,t) = 0, for t e dt„. The relation P(u,t) = 53 n°=1 Pn(w,t) must be fulfilled, and its physical meaning is obvious. TFD function can be discrete in time and frequency P(u^, tn)
Pk,n(d,t) — P (idk,tn ), idk — 1 < id < idk , tn — 1 < t < t>,
Pk,n(^,i) — 0, id < dk — 1, ^ > ^k , i < in — 1 , t > tn
k0 n0
P(d,t) = ^ ^ Pk,n(d,t)
k=1n=1
At the same time we introduce the normalized TFD functions P(uk,tn) and formulate the restricting set for the indices (k,n) e S and the optimization task
(k° ,n°) = arg{ max P (idk ,tn)}
k,n^S
4. Formulation of the Problem of Recognition of Anomalies in Time Series
We will use the standard methods of decision-making, which is usually applied to solving problems of statistical radio engineering [Levin, 1989]. Let us use the method, based on pattern recognition in frequency domain [Joswig, 1990; Savchenko, 1997]. Let us assume that a local reference TFD function P0(td, t) is independent of a local interval and is determined by a vector of parameters P0. We introduce the functional p0(Pn (d,t), P0(id,t)) = p0(Pn (id,t),P0), d e Q, t e dtn, equivalent of the distance between the evaluation of a local TFD function and local reference TFD function, which is accepted as a critical functional.
In this case the d.m.p. is based on the calculation of decisive functional for verification of the following conditions: if p0(Pn ,P0) = 1 - then a local interval dtn belongs to an anomaly, if p0(Pn, P0) = 0 - a local interval dtn is located within a non-anomalous section of a record.
Evaluation of boundary coordinates of anomalies is based on the analysis of the d.m.p. results for a sequence of local intervals and is probabilistic in nature. The decision whether a local interval belongs to an anomalous section or not, in practice, is always associated with errors: sometimes it is decided that an interval belongs to an anomalous section, and in reality it doesn’t, and vice versa.
Formulation of the recognition problem comprises two components.
• Evaluation of coordinates of the boundaries of anomalies
P(dk ,tn) — P (dk ,tn)/P (dk° ,tno )
Let us consider the signal x(t), for which the definition of anomalous and non-anomalous time series will be based on local TFD functions’ evaluations P° (d,t). We assume that anomalous and non-anomalous time series are determined by evaluations of local TFD functions with different spectral characteristics. In practice, these situations are fairly common.
Let the local intervals dtn1, dt„°, dtn3 correspond to the evaluations of local TFD functions P°1 — P°1(id,t), Pn2 = Pn2(u,t), Pn3 — Pn3(^,t). We suppose that the interval dtn1 is an anomaly and intervals dt„°, dtn3 are non-anomalous. We introduce, in the conventional sense, the distance functional p(•, •) between the evaluations of TFD functions that belong to different time series of records. We assume that for local intervals, belonging to non-anomalous series, TFD functions’ evaluations are almost identical; this means that p(P%°, Pn3) « 0. For local intervals, belonging to anomalous and non-anomalous series, the TFD functions differ significantly; this means that p(P„1,Pn°) > 0,
p(Pn1,Pn3) > 0.
Thus, the solution to the problem of recognition of anomalies in HPS signals is based on a calculation of local TFD and distance functional. We use the differences in the spectral parameters of signals in anomalous and non-anomalous intervals of time series.
(t°u,t°2l), I — 1, 2,... ,L°
• Evaluation of a probability of anomalies in the intervals (t°u,t°2l)
Pi , I — 1, 2,... , L°
5. Evaluation of TFD Functions
Let us proceed from a continuous to a discrete argument. Let us represent a time series of discrete signal as
y(Ti) — y(i), No < i < Nf
where N0, Nf are the numbers of the initial and final points of an initial signal x(t), T - a quantization step. We introduce a system of local intervals. Let dN be a number of points in a local interval and the number of the local boundary points of intervals satisfy the equations
Nln — No + (n — 1)dN
N°n — N1n + dN — 1, n — 1,...,no (1)
5 " , 10
0 0
Figure 6. TFD function for LR-disturbance.
where n0 - number of local intervals, n0 = ent((N/ — N0)/dN), function ent(-) - the integer part. Points, satisfying the condition N1n < i < N2n belong to the local interval n.
Let us consider the computing of sliding discrete Fourier transforms (DFT) [Lyons, 2007] for an initial signal. Each local interval will correspond to a N-points DFT-sliding interval. For the DFT-sliding interval with the step of sliding dN we assign the boundary points N1n, N2n. The interval corresponds to the local interval n
N1n = N1n — (N — dN )/2, N2n = N1n + N — 1 (2)
Let us formulate observations in the DFT-sliding intervals
y n( S) = y( N1n +S ), s = 0,... , N — 1, n =1,...,no
We write the expression for the coefficient of the DFT for the signal points in the n-th sliding DFT-interval and we apply the weighted Hanning window W(s),s = 0,... , N — 1 [Getmanov, 2010]
1 n—1 2ty( k — 1)s
°n(k) = N ^ yn(S)W(s)e ' N
s = 0
k =1,...,N/2, n = 1,...,n0 (3)
The index k defines the discrete DFT-frequencies Aw =
2-k/NT, Uk = Auc(k — 1), k = 1,... ,N/2. We take the expression
P°(cok, tn) = C* (k)Cn(k), k = 1,...,N/2
n = 1,...,no, P°(k,n) = P°(idk ,tn) (4)
as the evaluation P° (k, n) of a local TFD function in points k,n, corresponding to discrete frequencies Uk and local intervals tn.
Figure 7. TFD function for TW-disturbance.
6. TFD Functions Evaluation for HPS Signals with LR and TW-Disturbances
Let us analyze a HPS signal and consider a kind of TFD functions for sections with LR and TW-disturbances.
Figure 6 shows the evaluation of the normalized TFD function P°(k,n), corresponding to a case of LR. The interval 6,704,150 < i < 6,704,250 was taken where it was known a priori about a LR-disturbance; we used the N — 32-dimensional DFT for calculating the TFD with sliding step dN — 2. Figure 5 shows that the values of evaluated TFD are concentrated in a high-frequency area for indices k « 10 — 16.
Figure 7 shows the normalized TFD function P°(k,n), corresponding to a TW case. The interval 6, 707,125 < i <
6, 707, 205 was examined, the TFD dimension N — 32, sliding step dN — 2. We see that the TFD values are concentrated in a low-frequency area for indices k « 1 — 5.
It can be concluded from Figure 6 and Figure 7 that LR-disturbances correspond to low-frequency TFD. The conclusions are in compliance with the physics of LR and TW-processes. This conclusion can apply, after analyzing other intervals, to the entire signal 6,000,000 < i <
7, 600, 000.
7. Decision-Making Procedures and Indication Functions
7.1. Let us consider a case of high-frequency TFD functions, corresponding to LR-disturbances. Let P1r, P2r, k0R, N, P3r be the configuration parameters for the proposed d.m.p. The first four parameters P1r, P2r, k0R, N determine a reference normalized TFD function, the parameter
P3r defines the maximal value for the TFD function evaluation. We calculate relations for n — 1,... ,n0, in order to formulate a correct d.m.p., k1R — k0R + 1
P1MR,m(n) — max P!°n (k,n)
1<k<k0 R
P2MR,m (n) = max P^ (k,n)
k1R<k<N/2
PsMR,m(n) = max P^ (k,n) (5)
1<k<N/2
The proposed d.m.p., providing that the n local interval belongs to an anomalous section with a LR-disturbance, is related to computing the logical expression Pr(ti) using (5) for n = 1,... ,n0
Pr(ti) = (P1MR,m(n) < Pir) |^| (p2MR,m(n) >
P2R) |^| (PriMR,m(n) > P'3r) (6)
The formulae (6) have a quite clear physical meaning In this case the d.m.p. is to check the joint implementation of the inequalities (6) and calculation of the indicator coefficient Ir(ti): if the logical expression Pr(ti) (6) is correct, then the value Ir(ti) = 1 is accepted; otherwise - Ir(ti) = 0.
The sequence Ir(n), n = 1,...,n0, consisting of zeros and ones, is divided into rn0 indication intervals of equal size, equal to (INr, rn0dNR = n0. Initial and final points of an indication interval with the number rn can be calculated in the obvious way:
N1 R,(m) = N0 + dNR(m — 1)
N2R.(m) = N1R + dNp, m = 1, 2,...,mo (7)
Using (7) and the values Ir (n) the LR-indication function a0R(m) is calculated for each indication interval
1 N2R(m)
aoR.(m) = —— ^2 Ir(h), m = 1,2,...,mo (8)
R n = NlR(m)
7.2. Similarly, let us consider a case of low-frequency TFD functions, corresponding to TW-disturbances. Let P1t, P2t, k0T, N, P3t be setup parameters. For n = 1, ...,n0 we calculate the relations for TW, k1T = k0T + 1
P1MT,m(n) = max Pm (k,n)
1<k<kQ t
P2MT,m(n)= max?m (k,n)
kiT <k<N/2
P3MT,m(n) = max Pm (k,n) (9)
1<k<N/2
The proposed d.m.p. related to the fact that the local interval n belongs to an anomalous section with a TW-disturbance is implemented on the basis of calculating the local expression PT(n) for n = 1,... ,n0
Pt (n) = (p1MT,m(n) > p1T) ^ (p2Mt,m(n) <
P2T) ^ (Pr3MT,m(n) > P3T) (10)
The formulae (10), same as (6), are quite obvious and can be used for computing the indication coefficient IT(n): if the logical expression PT(n) (10) is true, then IT (n) = 1; otherwise - IT (n) = 0. The resulting sequence IT (n), n =
1,... ,n0 is divided into dNT, m0dNT = n0 - dimensional m0 indicator intervals. Initial and final points of the indication interval rn can be calculated in a similar way (7):
N1T (rn) = N0 + dNT (rn — 1)
N2T (rn) = N1T + dNT, rn = 1, 2,...,rn0 (11)
Applying (7) and IT (n) we calculate the TW-indication
function a0T (rn) for each indication interval
1 N2T(m)
O.0T (m) = It (n), m = 1, 2,...,mo (12)
n=^i^ (m)
The indication functions a0R(rn), a0T(rn) can be interpreted as the probability of the presence of LR and TW-disturbances at the indication interval m. Calculation of indication functions is connected with calculation of indicator functions associated with probability errors - possible false detection and omissions (probabilistic errors).
TFD functions evaluation (paragraph 5), implementation of d.m.p. and calculation of indication functions (paragraph 7) is made on the basis of STARTS-1 algorithm.
8. Non-Linear Amplitude and Time-Latitude Filtering of Indication Functions
To reduce the size of the text, the index “p” in the formulae will take two values: 1. “p”= R; 2. “p”= T, and correspond to the cases of LR and TW-disturbances.
The initial data for the implementation of the nonlinear filtering algorithm - the d.m.p. results are presented as the indicator function a0p(m), rn e M, M = {rn : (1 < rn < rn0)} from (8) or (12). The purpose of nonlinear filtering of indicator functions is to improve the accuracy of the boundaries of anomalous intervals and probabilities evaluations.
The proposed algorithm for nonlinear filtering of indicator functions is based on threshold filtering algorithms and algorithms for filtering pulse sequences with pulse modulation, known from radio engineering [Baskakov, 2005; Gonorovsky, 2006]. To a certain extent, this algorithm is heuristic, it requires a selection of configuration parameters and possible improvements. The implementation of nonlinear filtering is based on the STARTS-2 algorithm.
8.1. The first stage of the algorithm consists of realization of the nonlinear amplitude filtering of indicator function and is related to removing the small values from the sequence a0p(m), rn e M, M = {rn : (1 < rn < rn0)}. Removal of small values is done on the basis of comparison with the threshold a0p, and the sequence’s formation a1p(rn), rn e M
o.1P(rn) =0, if aos(rn) < aop
a1p(rn) = a0p(rn), if a0$(rn) > a0p,rn e M (13)
The threshold’s value parameter is adjustable and its correct choice ensures the optimal correlation between the probabilities of detection and error.
8.2. The algorithm’s second stage relates to implementation of non-linear time-latitude filtering of indicator function and is related to removing sites with a small amount of zero and positive values in the sequence from the first stage
a.1P(m).
For implementation of the second stage we use a function ^{-, •}, which calculates bordering points of intervals with zero and positive values of sequence elements of a certain input sequence y(rn),rn e M and the corresponding numbers of these intervals. In general, we get
^{y(rn), m e M} = {(m1,iom2,io), 10 = 1, ...,Lo,
(rn1,i1,rn2,n), 11 = 1,...,Ln}
where (rn1,i0,rn2,i0) - bordering points of intervals with zero elements, L0 - number of intervals with zero elements, (rn1,i1,rn2,i1) - bordering points of intervals with positive elements, L1 - number of intervals with positive elements. The calculations on the function ^{-, •} are done with the help of the CABORD (Calculation of Borders) algorithm.
Let us apply the function ^{-, •} to a1p(rn), we find the bordering points of intervals with positive elements and the corresponding numbers of these intervals
fy{a1P(m),m e M} = {(m1,n,m2,n), 11 = 1,...,L,1} (14)
Further, on the basis of (14) the set L1 of the intervals’ numbers can be found, which satisfy the inequality (rn2,i1 — m1,i1) < drn1p, where drn1p is the setup parameter. Then the corresponding set M1 is formed on the basis of indices rn
L1 = {/1 : (rn,2,n — m1,n < drn1P, 11 = 1,. .. , L1)}
M1 = {rn : (m1,i1 < rn < m2,i1, 11 e L1)} (15)
on its basis the interim sequence a2p(rn) is formed, where single (or rare) positive values in a1p(rn) are replaced by zeros
a2p(rn) = 0, rn e M1;
a2p(rn) = a1p(rn), rn e M/M 1 (16)
After that intervals with a small number of zeros are removed from the sequence a2p(rn). For this purpose the function ^{-, •} is used again for a2p(rn), rn e M, for which the borders of intervals with zero elements are defined
^{_U2P(m),m e} = ((m1,io,m2,io), 10 = 1,. ..,Lo) (17)
and from (17) the corresponding sets L0, M0 are found, taking into account the setup parameter dm2p.
L0 = {/0 : (rn2,io — m,1,io < dm,2P, 10 = 1, ...,Lo)}
M0 = {m : (m1,i0 < m < m2,i0,l0 e ¿0)} (18)
The result of the non-linear filtering is represented as the sequence £p(rn)
(rn) = a0p(rn), rn e M0
£°(m) = a2p(rn), m e M/M0 (19)
9. System of Algorithms of Recognition of Anomalies
The solution of our problem of recognition of abnormal areas is based on algorithms for evaluation of the boundaries of abnormal areas and calculation of the probabilities.
The result of the filtering £,°(m) is used, which is transformed by the function ^{-, •}, with the purpose to find the sequence of borders of intervals with positive elements (m°,i,m°,i), I — 1,...,L°
^{£p(m),m € M} — (m°,i, rn°,i), I — 1,...,L°,
i — No + (m — 1)dNo, (n,i,i°,i), I — 1,...,L° (20)
Evaluation of probabilities of anomalies in the intervals (rn1i, rn°i), I — 1,... ,L°, is done by the following obvious formulae
m2,i
Pip — ( aoP(rn))/(rn°°:i — m°1,i + 1),
m=m'2 i
l — 1,...,L°, aip — pip (21)
Figure 8. Graph of indicator function aR(i) of recognition of intervals with LR-disturbances.
Calculation of the indicator function ap(i) for No < i < Nf is done by using (21)
aip(i) = piv for i e (*°,i,*°,i), aip(i) = 0 for i e (i°,i ,i°,i)J = 1,... ,L°,
aP(i) = Yl aip(i) (22)
i=i
Thus, the STAN method for the present task of recognition led to a sequential application of the STARTS-1 and STARTS-2 algorithms. Let us describe the procedure.
In accordance with the STARTS-1 algorithm local TFD functions are evaluated and implementation of the d.m.p. coming to actions 1-3, which are related to:
• Formation of systems of local and DFT - sliding intervals, formulae (1), (2);
• Evaluation of local TFD functions, formulae (3), (4);
• D.m.p. implementation and calculation of indicator functions, formulae (5)-(8), (9)-(12).
For the STARTS-1 algorithm we list: input variables -N0,Nf, dN, N, y(i),N0 < i < Nf; setup parameters -PiR, P2r, koR, P3r; output variables - aoR(to), aor(to), to = 1, 2,... , to0.
On the basis of the STARTS-2 algorithm the non-linear filtering of indicator functions is done, which is reduced to steps 1-3, which in their turn comprise:
• Amplitude filtering, formulae (13);
• Time-latitude filtering, formulae (14)-(19);
• Calculating the borders of evaluations of anomalies and probabilities’ evaluations, formulae (20)-(22).
For the STARTS-2 algorithm we list: input variables - a0p(rn), to = 1,2,.. . ,m0; setup parameters - a0p, dm, 1p, drn2p; output variables - £°(to), to = 1, 2,...,m0, (mli,to,2,i), I = 1,...,L°, Pip, I = 1,...,L°, ap(i). The STARTS-2 algorithm uses the results of the application of the CABORD algorithm.
10. Testing of System of Recognition Algorithms
We considered the filtered HPS signal y(i) = y(Ti) from NGDC at a fixed time interval; No = 600,001, Nf =
7, 600, 000 - numbers of the initial and final point of this interval. The signal was used for testing the developed system of the STARTS-1 and STARTS-2 sequential algorithms for recognition of anomalies in time series with LR and TW-disturbances. The results of the algorithms’ application were compared to the control data, received from NGDC.
10.1. Recognition of time series with LR-disturbances was implemented. A fixed time interval was divided into no = 3, 500, 000 dN = 2-dimensional local intervals. Each of the local interval corresponded to DFT-sliding intervals with dimension N = 32. dNR = 20 - number of points of an indicator interval, to0 = 175, 000 - number of indicator intervals. The fixing parameters for d.m.p. took the values: PiR = 0.10, P2 R = 0.3, P3 r = 0.0005, koR = 10. For the setup parameters of non-linear filtering we took the values aRo = 0.4, dmir = 2, dTO2r = 2.
Figure 8 shows the results of the algorithm’s application
Table 1. The Data on Recognition of LR-Disturbances
ii NDBC Î2 NDBC ii GC RAS Î2 GC RAS aiR
1 660,314 660,339 660,3216 660,382 0.77
2 998,466 998,824 998,511 998,562 0.97
3 2,235,851 2,235,915 2,235,871 2,235,912 0.87
4 2,394,185 2,394,337 2,394,221 2,394,262 0.98
5 3,055,796 3,057,107 3,056,341 3,05,6402 0.99
6 3,057,394 3,059,870 3,056,511 3,056,572 0.98
7 3,080,259 3,080,404 3,080,251 3,080,322 0.98
8 3,709,684 3,709,735 3,709,691 3,709,732 0.45
9 4,416,084 4,416,257 4,416,081 4,416,120 0.97
10 5,100,460 5,100,557 5,100,461 5,100,502 0.44
11 6,704,179 6,704,246 6,751,311 6,751,352 0.66
12 6,985,705 6,985,835 6,985,771 6,985,812 0.53
13 7,571,284 7,571,402 7,571,351 7,571,402 0.89
in the form of a graph of the indicator function ur (i), representing a set of vertical lines of different height.
Table 1 comprises the data on recognition of LR-disturbances. Columns 2, 3 show the coordinates of the bordering points of intervals with LR-disturbances from NGDC (i1NGDC, ¿2NGDC), which were calculated with the help of a semi-automatic NGDC-algorithm, based on the analysis of signals in time domain. Columns 4, 5 show the coordinates of the bordering points of intervals with LR-disturbances, calculated by the system of algorithms, proposed by Organization of the Russian Academy of Sciences Geophysical Center (GC RAS). Column 6 shows the evaluations of probabilities of LR-disturbances aiR in the intervals with calculated bordering points (i1li,i2,i), I = 1,..., Lq, Lq = 13.
We note that the intervals #5 — 7 with LR-disturbances, according to Table 1, have been recognized (their time resolution was realized); however, in Figure 7 these intervals are shown together. The proposed algorithm has implemented a sufficiently precise recognition of almost all the marked LR-disturbances. It should be noted that in interval #11, in the points 6, 704,179 < i < 6, 704, 246, mentioned by NGDC, there was really a LR-disturbance, which was weaker than another LR-disturbance, located in the interval
6, 751, 311 < i < 6, 751, 352.
10.2. Recognition of time series with TW-disturbances. A fixed time interval was divided into nQ = 1, 750, 000 local intervals of the dimension dN = 4. Each of the local intervals corresponded to DFT-sliding intervals of the dimension
Figure 9. Graph of the indicator function ar (i) of recognition of intervals with TW-disturbances.
Table 2. The Data on Recognition of TW-Disturbances
ii NDBC i2 NDBC ii GC RAS Î2 GC RAS am
1 661,933 665,749 662,000 664,782 0.91
2 1,000,907 1,002,656 1,000,500 1,001,032 0.51
3 6,706,714 6,708,195 6,706,750 6,707,532 0.61
4 7,574,954 7,575,038 7,577,250 7,577,282 0.99
N = 32. We assumed that dNr = 25 was the number of points of the indicator interval, mQ = 70, 000 - the number of indicator intervals. The setup parameters for d.m.p. took the values: P-\_t = 0.5, P2t = 0.125, P3t = 0.004, kQT = 5. The values of setup parameters of non-linear filtering were a-TQ = 0.4, dm1T = 2, dm2T = 2.
Figure 9 shows the graph of the indicator coefficient ar (i) as a set of vertical lines of different height. Their intersection with the horizontal axis of coordinates determines the time boundary points of intervals with TW-disturbances.
Table 2 shows the data on recognition of TW-disturbances. Columns 2, 3 contain the coordinates of the bordering points of TW-intervals, determined on the basis of the semi-automatic recognition algorithm, based on the analysis of signals of time series. The algorithm was developed by NDBC and the data was transferred to GC RAS. Columns 4, 5 show the coordinates of the bordering points of TW-intervals, obtained on the basis of the algorithms of recognition, developed by GC RAS. Column 6 contains the evaluations of probabilities of recognition of TW-disturbances uti in the intervals with calculated bordering points (i1,i, *2,0, I = 1,... , Lq, Lq = 4.
11. Conclusion
The proposed STAN method for recognition of time series with LR and TW-disturbances, if we analyze the results of Figure 8 and Figure 9 and Table 1 and Table 2, works quite well. Thus, for 13 + 4 = 17 disturbances the proposed algorithms allowed only one false disclosure and one omission, that allows us to conclude, at a first approximation, that the level of false disclosures and omissions was equal to ~ 5.8%.
The effectiveness of the STAN method of recognition of time series with LR and TW-disturbances in terms of minimizing the probability of false detections and omissions, as well as possible refinements of the boundary points of anomalies can be improved by optimizing the choice of setup parameters. Due the adequate mathematical apparatus of TFD functions, the proposed STAN method and the corresponding algorithms have the potential to provide detection of weak LR and TW-disturbances.
The proposed STAN method of recognition can be used to automate the viewing of archives records of observations from the BSS network and to ensure the efficient solution of the tsunami warning problem.
The developed method of STAN-recognition is quite versatile and enables the solution of problems of recognition
of anomalies in time series of geophysical data of different nature.
Acknowledgments. The author considers it his pleasant duty to express his gratitude to Academician A. D. Gvishiani for helpful discussions of the results of this work and to the scientists of the National Geophysical Data Center of the U.S. K. Stroker, G. Moungoff for providing very useful informational support.
References
Agayan, S. M., S. R. Bogoutdinov, A. D. Gvishiani, E. M. Graeva, J. Zlotnicki, M. V. Rodkin (2005), Study of signal’s morphology on the basis of fuzzy logic algorithms, Geophysical Studies, 1,
1, 143-155.
Baranov, S. V. (2007), Application of wavelet-transform for automatic detecting of seismic signals, Fizika Zemli, 2, 83-94.
Baskakov, S. I. (2005), Radio Engineering Networks and Signals, Moscow, Vyshaya Shkola, 462.
Bashilov, I. P., et al. (2008), Ocean-bottom geophysical observatories — methods of construction and fields of application, Nauchnoe Priborostroenie, 18, 2, 86—97.
Boashash, B. (2003), Time Frequency Signal analysis and processing, Elsevier Science, 770.
Bogoutdinov, S. R., A. D. Gvishiani, S. M. Agayan, A. A. Soloviev, E. Kihn (2010), Recognition of disturbances with given morphology in time series I. spikes on INTERMAGNET magnetic records, Fizika Zemli, 11, 99—112.
Chebrov, D. V., A. A. Gusev (2010), Automatic definition of parameters of tsunami generating earthquakes in the russain Far East in real time: Algorithms and software, Seismicheskie Pri-bory, 46, 3, 5—21.
Cohen, L. (1989), Time-frequency distributions: review, TIIER, 77, 10, 72—120.
Dykhan, B. D., V. M. Zhak, E. A. Kulikov (1981), First tsunami registering in open ocean, DAN, 257, 5, 1088—1092.
Dzienovski, A., S. Bloch, M. Landisman (1969), Technique for the analysis of transient seismic signals, Bull. Seismol. Soc. Amer., 59, 427—444.
Getmanov, V. G. (2010), Digital Processing of Signalds, Moscow, NIYaU MIFI Publishers, 248.
Getmanov, V. G., A. D. Gvishiani, K. Stroker (2011), Recognition of P-waves disturbances in observations of specialized bottom stations. Deep structure, geodynamics, interpretation of geophysical fields, Shestye Nauchnue Chtenija Pamyati Y. P. Bulashevicha, Materialy Konferencii, Ekaterinburg, URO RAN, 76—79.
Gonorovsky, I. S. (2006), Radio Engineering Networks and Signals, Textbook, Moscow, Drofa, 719.
Gvishiani, A., J. Dubois (2002), Artificial Intelligence and Dynamic Systems for Geophysical Applications, Paris, SpringerVerlag, 350.
Gvishiani, A. D., S. M. Agayan, S. R. Bogoutdinov (2008), Definition of anomalies in time series by methods of fuzzy recognition,
Doklady RAN, 421, 1, 101—105.
Hlawatch, F., F. Auger (2008), Time-Frequency Analysis: Concepts and Methods, London, ISTE and Wiley, 440.
Joswig, M. (1990), Pattern recognition for earthquake detection, Bulletin of the Seismological Society of Placecountry-Region America, 80, 1, 170-186.
Kedrov, E. O., O. K. Kedrov (2006), Spectral-time method of identification of seismic phenomena at distances 15-40°, Fizika Zemli, 5, 47-64.
Kedrov, O. K., L. A. Polikarpova, G. M. Steblov (1998), Algorithm of recognition of weak short-period seismic signals on the basis of time-frequency analysis of three-component records in real time regime, Izv. RAN, Fizika Zemli, 2, 30-45.
Kulikov, E. A., F. Gonzales (1995), Restoration of tsunami waveforms at the source by measuring vibration level of ocean by remote hydrostatic pressure sensor, DAN, 344, 6, 814-818.
Lyons, R. (2007), Digital Processing of Signals (Translated from English), Moscow, Binom-Press, 656.
Lander, A. V., A. L. Levshin, V. F. Pisarenko (1973), Spectraltime analysis of oscillations, Vychislitelnaya Seismologija, 6, 3-27.
Levin, B. R. (1989), Theoretical Basics of Statistical Radio Engineering, Moscow, Radio I Svjaz, 665.
Levin, B. R., M. A. Nosov (2005), Physics of Tsunami and Related Phenomena in Ocean, Moscow, Janus-K, 360.
Meining, C., S. Stalin, A. I. Nacamura, F. Gonzalez (2005), Technology development in real-time tsunami measuring, monitoring and forecasting, Proc. IEEE Oceans Conf., Sept. 2005, 2, 1673—1679.
Poplavsky, A. A., E. A. Kulikov, L. N. Poplavsky (1988), Methods and Algorithms of Operative Tsunami Forecasting, Moscow, Nauka. 128.
Savchenko, V. V. (1997), Recognition of random signals in frequency domain, Radiotekhnika I Elektronika, 42, 4, 426—429.
Tompson, S. R. (2009), Sound Propagation Considerations for a Deep-Ocean Acoustic Network, US Naval Academy, 63.
V. G. Getmanov, Geophysical Center RAS, 3 Molodezhnaya Str., 119296 Moscow, Russia. ([email protected])
