CC BY
10
UDC 004.891:550.3 DOI: 10.25513/2222-8772.2018.3.91-100
SOFT COMPUTING IDEAS CAN HELP EARTHQUAKE
GEOPHYSICS
Solymar Ayala1
Doctoral Student, e-mail: [email protected] Aaron Velasco1 Ph.D. (Geosciences), Professor, e-mail: [email protected] Vladik Kreinovich2 Ph.D. (Phys.-Math.), Professor, e-mail: [email protected]
1 Department of Geological Sciences 2Department of Computer Science
University of Texas at El Paso, Texas 79968, USA
Abstract. Earthquakes can be devastating, thus it is important to gain a good understanding of the corresponding geophysical processes. One of the challenges in geophysics is that we cannot directly measure the corresponding deep-earth quantities, we have to rely on expert knowledge, knowledge which often comes in terms of imprecise ("fuzzy") words from natural language. To formalize this knowledge, it is reasonable to use techniques that were specifically designed for such a formalization — namely, fuzzy techniques. In this paper, we formulate the problem of optimally representing such knowledge. By solving the corresponding optimization problem, we conclude that the optimal representation involves using piecewise-constant functions. For geophysics applications, this means that we need to go beyond tectonic plates to explicitly consider parts of the plates that move during the earthquake. We argue that such an analysis will lead to a better understanding of earthquake-related geophysics.
Keywords: earthquake geophysics, soft computing, Haar wavelets.
1. Specifics of Data Processing in Earthquake Analysis (and in Geophysics in General)
Earthquake analysis is important. Earthquakes can be devastating. It is therefore important to gain as much understanding about the corresponding geophysical processes as possible; see, e.g., [1,7].
Usual approach to data processing. A good understanding means that we know, for each location (x,y) and for each depth z, what is the density p at this location and this depth, what are the mechanical properties of the material at this 3-D location u = (x,y,z), what are the stresses at this 3-D location. In other words, we need to find the corresponding functions like p(u) = p(x,y,z).
How can we describe a function in a computer? A usual way is to select a basis of functions
ei(u),e2(u),..., (1)
so that each desired function f (u) can be represented as a linear combination of the basis functions
f (u) = ci • e\(u) + C2 • u2(u) + ..., and then represent the desired function f (x) by the corresponding coefficients
(Ci ,<*,...). (2)
In principle, we can consider different bases, but it is usually convenient to orthonormalize them, i.e., to consider linear combinations
i
CM = ^ Cij • e3 (u)
3 = 1
for which, for all i and j, we have
J (e™(u))2 du = 1
and
J e™(u) • ef(u) du = 0 when i = j. In this case, the desired coefficients Ci can be obtained by using a simple formula
Ci = j f (u) • e™(u) du. (3)
Thus, without losing generality, we can safely assume that the basis (1) is orthonormal.
The most widely used examples of such bases are:
sines and cosines, and • wavelets; see, e.g., [2,4,6,11]. For sines and cosines, the expansion into the corresponding basis is known as Fourier transform. For wavelets, the transformation from the original function f (u) to the coefficients q is known as the wavelet transform.
It is important to select an appropriate basis. It is known that selecting an appropriate basis can drastically improve the quality of the data processing results.
For example, in many cases, wavelet analysis has led to interesting discoveries that were not possible when Fourier analysis was used to process the corresponding data.
It is therefore very important, in each practical situations, to select the most appropriate basis.
What we plan to do in this paper: main idea. In this paper, we provide arguments for selecting the most appropriate basis for earthquake-related analysis. In this analysis, we use the specific features of the geophysical data processing.
Specifics of geophysical data processing. In comparison with most other data processing situations, geophysical analysis has two important specifics.
First, in most data processing situations, we have continuous functions. For example, when we control a vehicle, its location continuously depends on time. In contrast, in geophysics, there are clear discontinuities: as we go deeper,
we have an abrupt transition between different layers. The second difference is that in most other data processing situations, we can determine the ground truth, i.e., the actual values of the corresponding quantities. In geophysics, our ability to get the ground truth is very limited: up to a certain depth, we can drill a borehole and find out what are the actual properties, but at larger depths, this is not practically possible.
Why soft computing. Since we cannot determine the actual values to check different models, we have to rely on expert knowledge to decide which model works better.
Expert knowledge rarely comes in precise terms, it usually comes in terms of imprecise ("fuzzy") words. To describe the corresponding knowledge in precise terms, it is therefore reasonable to use techniques specifically designed to handle such knowledge — namely, the techniques of fuzzy logic; see, e.g., [3,5,8-10,12, 13].
2. Analysis of the Problem
Main idea. Since the values f (u) comes from expert estimates, they come with a fuzzy uncertainty. In other words, for every u, we have fuzzy information about the difference Af (u) =f f(u) — f (u) between:
• the expert estimate f(u) and
• the actual (unknown) value f (u). In precise terms, this means that:
• we do not know the probabilities of different possible values of Af (u), but
• we have a membership function p(Af) that describes, for each possible value Af, the degree to which this value is possible.
Since we have no reason to assume that the estimation errors are positive or negative, it is reasonable to assume that the degree of possibility of each value Af does not depend on its sign: p(—Af) = p(Af).
Once we have selected the basis ei(u), we will then transform the estimate for f (u) into the sequence of the corresponding coefficients q.
• Since the values f (u) are known with uncertainty,
• as a result, we can only determine the coefficients q — and thus, the corresponding term d • ei(u) — with uncertainty.
A reasonable idea is to select the basis eH(u) for which this uncertainty in the term d • d(u) is the smallest possible.
Let us describe this idea in precise terms. When we process the expert estimates f(u), we get the following estimates Z for the coefficients cc
Z = J f(u) • ei(u) du. (4)
The actual (unknown) value q of the corresponding coefficient can be obtained if we apply the same procedure to the actual (unknown) function f (u):
Ci = J f (u) • ei(u) du. (5)
If we subtract (5) from (4) and take into account that the integral of the difference is equal to the difference of the integrals, we get the following formula for the inaccuracy Aci = Z — q:
Ac, = J Af (u) • ez(u)du. (6)
The inaccuracy in the product Ci • ei('u) is equal to the product Aci • ei('u). This value depends on the location u:
• for some locations u, the value lei(u)l is larger, so the inaccuracy is larger;
• for other locations u, the value |ei(u)l is smaller, so the inaccuracy is smaller. It is reasonable to minimize the worst-case inaccuracy
A Ci • max | ei(u)l = max | ei (u)l • A f(u) • ei(u)du. (7)
u u J
Here, each value A f(u) is a fuzzy number, so Ad is also a fuzzy number. In fuzzy logic, this fuzzy number is determined by the Zadeh's extension principle.
It is known that in general, computing the result Y = f(Xi,... ,Xn) of applying a function f(xi,..., xn) to n fuzzy numbers X-]^,... ,Xn can be described as follows:
for each a e (0,1], the a-cut Y(a) d=f {y : ^(y) ^ a] is equal to the range of the function f(xi,... ,xn) when each Xi takes values from the corresponding a-cut
Xi(a) d=f {xi : ^i(xi) ^ a]. In precise terms, we have
Y (a) = {f (xi,... ,xn) : xi e X- (a),... ,xn e Xn(a)].
Let us apply this general result to our formula (7). Since the membership function ^(Af) does not change if we change the sign of the difference Af, for each a, the corresponding a-cut is a symmetric interval. Let us denote this interval by
[-A(a), A(a)].
The expression (7) is a linear combination of all the values A f(u):
• when ei(u) > 0, this function is increasing in A f (u);
• when ei(u) < 0, this function is decreasing in Af (u).
Thus, when each value Af(u) takes all possible values from the interval [—A(a), A(a)], the largest possible value of the expression (6) is attained when:
• for those u for which d(u) > 0, the value Af (u) is the largest possible, i.e., Af(u) = A(a), and
• for those u for which ei(u) < 0, the value Af (u) is the smallest possible, i.e., Af (u) = -A(a).
In both cases, the largest possible value of the product Af (u) • d(u) is equal to A(a) • lei(u)l Thus, the largest possible value of the integral (6) is equal to
J A(a) • |ei(«)| du = A(a) • j |ej(«)| du.
Hence, the largest possible value of the integral (7) is equal to
A(a) • max lei(u)l • l^i(u)l du. u J
Similarly, we can show that the smallest possible value of the expression (7) is equal to
-A(a) • max lei(u)l • l^i(u)l du.
u J
Thus, Y (a) = [-y(a),y(a)], where
y(a) d=f A(a) • max lei(u)l • I^î(u)I du.
u J
The estimate for q is the most accurate when this interval is the narrowest possible, i.e., when the value
A(a) • max lei(u)l • l^i(u)l du u J
is the smallest possible.
In this product, the factor A(a) is given. So, the smallest possible value of the above product is attained when the product
max lei(u)l • l^i(u)l du (9)
u J
attains its smallest possible value. Hence, we arrive at the following optimization problem.
Resulting optimization problem. Among all possible functions ei(u) for which f ef(u) du = 1, we need to find the function with the smallest possible value of the product (9).
Analysis of the resulting optimization problem. We always have el(u) = lei(u)l2. Thus, / lei(u)l2 du = f e?(u) du = 1.
Also, for every u, we have ei(u) ^ maxv |ei(v)l Hence,
IZî(U)I2 ^ (max |ei(v)\^ • |ei(u)\
(10)
Integrating both parts of this inequality, we conclude that
1
J \ei(u)\2 du ^ max \ei(v)\ • J \ei(u)\ du
(11)
Thus, the product (9) that we want to minimize cannot be smaller than 1. One can easily check that when |ei(u)l = const for all u from the given region, we get exact equality in the formula (1) and thus, in formula (11).
So, when the absolute value |ei(u)l is constant, we attain the smallest possible value 1 of the desired product (9).
Vice versa, if at least for one value u, we have strict inequality in (10), we will have strict inequality in (11) as well. So, to attain the smallest possible value of the product (9), we must always have equality in the formula (1), i.e., we must always have the following equality:
When d(u) = 0, we can divide both sides of this equality by |e i(u)l and conclude that |ei(u)l = maxv |ei(v)l In other words, for every u:
• we either have ei(u) = 0
• or we have |ei(u)l equal to the largest possible value m = max |ei(v)l
V
So, we arrive at the following conclusion:
Resulting solution. Each function ei(u) from the geophysically optimal basis take only three values: 0, m, and —m, for some real number m> 0.
This means, in particular, that all the optimal basis functions are piecewise-constant.
Comment. Let us consider the geophysical meaning of this result.
3. Geophysical Meaning of Our Result
What does our result means in terms of earthquake analysis. An earthquake leads to a spatial shift at different locations. For catastrophic earthquakes, this shift can be in meters; for smaller earthquakes, we can have centimeters-size shift.
In general, we have a shift s(x, y) as a function of 2-D spatial coordinates. Our optimization result shows that the optimal way to analyze the empirical data about this shift is to represent it as a linear combination d• ei(u) of piece-wise constant functions ei(u). Such a linear combination is also piece-wise constant. Thus, what we need to do is to divide the whole area into several zones, in each of which the shift is fixed.
\ei(u)\2 = (max \e¿(w)\) • \ei(u)\
(12)
In geometric terms, this means that instead of considering each spatial location (x,y) by itself, we divide the whole region into parts, each of which moves as a whole (i.e., as a solid body).
How do we transform the observed shifts into this piece-wise constant presentation: an algorithm. When a function is piece-wise constant, it means that it attains finitely many different values. Let us sort these values into an increasing order: s1 < s2 < ... < sm.
Suppose at first that these values are given. In this case, we want to approximate the original function f (u) by a piece-wise constant function a(u) that takes values Si. For each u, the value a(u) is equal to one of the values s1,... ,sm. Thus, describing the function a(u) is equivalent to describing, for each i from 1 to m, the set Si of all the locations u for which a(u) = Si. These m sets should form a partition of the original domain S, i.e., we should have S1 U ... U Sm = S and Si n Sj = 0 for all i = j.
A natural idea is to use the Least Squares approach, i.e., to find such a function a(u) for which the integral /(f (u) — a(u))2 du attains the smallest possible value. One can easily check that the integral attains the smallest possible value if and only if for each u, we select the value a(u) e {s1,...,, sm} for which the value (f (u) — a(u))2 is the smallest possible. In other words, for each location u, as a(u), we take the value Si which is the closest to the original value f (u). In other words:
• we select a(u) = s1 if f (u) ^ Si + Si;
I + ( \ -f S1 + S2 / t( N ^ S2 + S3
• we select a(u) = s2 if —-— ^ j(u) ^ —-—; ...
for each i = 2,... ,m — 1, we select a(u) = Si if S%~1 —— ^ f(u) ^ S% ■
2 ^ 7 2
finally, we select a(u) = sm if Sm 1 — Sm ^ f (u).
We can repeat this procedure for different tuples s = (s1,..., sm). For each such tuple, we find the resulting mean square error
J min(f (u) — Si)2 du.
We then select a tuple s for which this mean square error attains the smallest possible value.
This immediately brings to mind tectonic plates. The above piece-wise description bring to mind the geophysical idea that the earth's surface consists of tectonic plates, solid bodies that move in relation to each other.
So, in the first approximation, our mathematical result leads to the very well-known plate tectonics idea.
Our result goes beyond plate tectonics. In the first approximation, our result simply leads to a well-known idea of plate tectonics. In this first approximation, the whole plate moves as a whole, the shift is exactly the same on all locations from this plate.
In practice, the shift is somewhat different in different locations on the same tectonic plate. To capture this difference and thus, provide a more accurate description of the corresponding geophysics, we therefore need to divide each affected plate into two (or more) different parts, with different shifts in each part.
This idea has geophysical sense. It is known that the major earthquakes are caused by the interaction of tectonic plates — that move relative to each other. As a result, all major earthquakes — and the vast majority of smaller earthquakes -happen at the boundaries between tectonic plates. Specifically, they happen at the convergent boundaries, where the plates move towards each other, accumulating a stress. This stress is released by an earthquake.
The above description is a first crude approximation to the corresponding geophysics, in which we can consider the whole plate as a solid body, in which all parts move the same way. In reality, different parts of the plate may accumulate the stress differently and move differently. As a result, some earthquakes only involve a part of the boundary between the plates. Depending on the size of this part, we can get earthquakes of different magnitudes.
Beyond piece-wise constant functions: geophysics-motivated idea. Solid bodies do not just shift, they can also rotate. So, a natural idea is to consider not only shifts, but also rotations of the parts of the plate. In this case:
• instead of approximating the measured values f(u) by a piece-wise constant function,
we approximate it by a piece-wise linear functions corresponding to shifts and rotations of different parts of each tectonic plate.
This can help in earthquake studies. In view of the above, to get a better understanding of the earthquake geophysics, it is important to analyze which parts of the plate are involved in different earthquakes, which parts have accumulated more stress and in which part, the stress has been released.
This idea is challenging. From the computational viewpoint, our idea is very challenging:
while we can relatively easily identify the boundary between the plates, where the big motion occurs,
it is much more challenging to identify the parts of the plate that are involved in an earthquake.
The reason why this identification is not easy is because we are interested in geophysical processes far away from the boundaries, where the earthquake-related motion is much smaller in amplitude and thus, much more difficult to detect.
We all need to work together to overcome these challenges. As of now, what we have are ideas and models.
Our preliminary results show that these ideas are promising, and we will continue working on them.
However, we think that it will be beneficial to publicize these ideas so that others can implement them, use them, improve them if needed — and thus, help
to get a better understanding of earthquake-related geophysics. Acknowledgments
This work was supported in part by the CONACYT Fellowship for Solymar Ayala Cortez and by the US National Science Foundation grant HRD-1242122.
The authors are greatly thankful to Dr. Throne Lay from the Earth and Planetary Sciences Department of the University of California at Santa Cruz and to Dr. Steven Harder from the Department of Geological Sciences of the University of Texas at El Paso for valuable discussions and encouragement.
References
1. Aki K., Richard P. Quantitative Seismology, Theory and Methods. University Science Books, Sansalito, California, 2009.
2. Addison P.S. The Illustrated Wavelet Transform Handbook: Introductory Theory and Applications in Science, Engineering, Medicine and Finance. CRC Press, Boac Raton, Florida, 2017.
3. Belohlavek R., Dauben J.W., Klir G.J. Fuzzy Logic and Mathematics: A Historical Perspective. Oxford University Press, New York, 2017.
4. Chandrasekhar E., Dimri V.P., Gadre V.M. (eds.) Wavelets and Fractals in Earth System Sciences. CRC Press, Boca Raton, Florida, 2014.
5. Hajan A., Styles P. Application of Soft Computing and Intelligent Methods in Geophysics. Springer, Cham, Switerland, 2018.
6. Han B. Framelets and Wavelets: Algorithms, Analysis, and Applications. Springer, Cham, Switzerland, 2017.
7. Havskov J., Ottemoller L. Routine Data Processing in Earthquake Seismology: With Sample Data, Exercises and Software. Springer, Dordrecht, Heidelberg, London, New York, 2010.
8. Klir G., Yuan B. Fuzzy Sets and Fuzzy Logic. Prentice Hall, Upper Saddle River, New Jersey, 1995.
9. Mendel J.M. Uncertain Rule-Based Fuzzy Systems: Introduction and New Directions. Springer, Cham, Switzerland, 2017.
10. Nguyen H.T., Walker E.A. A First Course in Fuzzy Logic. Chapman and Hall/CRC, Boca Raton, Florida, 2006.
11. Nicholas P. Framelets and Wavelets: Algorithms, Analysis, and Applications. Cambridge University Press, Cambridge, UK, 2017.
12. Novak V., Perfilieva I., Mockor J. Mathematical Principles of Fuzzy Logic. Kluwer, Boston, Dordrecht, 1999.
13. Zadeh L.A. Fuzzy sets // Information and Control. 1965. V. 8. P. 338-353.
ИДЕИ МЯГКИХ ВЫЧИСЛЕНИЙ МОГУТ ПОМОЧЬ ГЕОФИЗИКЕ
ЗЕМЛЕТРЯСЕНИЯ
Солимар Айяла1
аспирант, e-mail: [email protected]
Аарон Веласко1 к.г.-м.н., профессор, e-mail: [email protected] В. Крейнович2
к.ф.-м.н., профессор, e-mail: [email protected]
1 факультет геологических наук 2факультет компьютерных наук Техасский университет в Эль Пасо, США
Аннотация. Землетрясения могут быть разрушительными, поэтому важно получить хорошее представление о соответствующих геофизических процессах. Одна из проблем в геофизике заключается в том, что мы не можем непосредственно измерять соответствующие глубинные величины, мы должны полагаться на экспертные знания, знания, которые часто выражены неточными («нечёткими») словами естественного языка. Чтобы формализовать эти знания, разумно использовать методы, специально предназначенные для такой формализации, а именно, нечёткие методы. В этой статье мы формулируем задачу оптимального представления такого знания. Решая соответствующую задачу оптимизации, мы заключаем, что оптимальное представление включает использование кусочно-постоянных функций. Для приложений геофизики это означает, что нам нужно выйти за пределы тектонических плит, чтобы подробно рассмотреть части плит, которые движутся во время землетрясения. Мы утверждаем, что такой анализ приведёт к лучшему пониманию геофизики, связанной с землетрясением.
Ключевые слова: сейсмическая геофизика, мягкие вычисления, вейвлеты Хаара.
Дата поступления в редакцию: 30.06.2018
