RUSSIAN JOURNAL OF EARTH SCIENCES VOL. 9, ES1001, doi:10.2205/2007ES000217, 2007
Change of epochs in Earth sciencesV. N. StrakhovInstitute of Physics of the Earth, Russian Academy of Science, Moscow, Russia ContentsAbstract[1] The present work tells about the author’s views regarding the two subjects. The first subject — the emergence of the new epochs in the Earth sciences, namely the epochs of geophysics and geoinformatics (previously there were two epochs in the Earth sciences — the epoch of geography and the epoch of geology). The second subject — the formation of the principally new theory of potential fields interpretation (gravitational and magnetic anomalies). Introduction[2] Only in 2004, reflecting on the future of geophysics and the essence of the future geoinformation science (geoinformatics), I realized that two epochs took place in the past evolution of Earth sciences, and the third epoch started developing since the end of the 20th century (tentatively, since 1995). [3] The first epoch of the Earth science evolution is the epoch of exploring the Earth's face (i.e. the structure of its surface). This epoch, which started from famous geographic discoveries and continued until the mid-19th century, should be called a geography epoch. It is in this epoch that geographic and geodetic sciences developed and geoinformation started it to be mapped cartographically. [4] The second epoch of the Earth science evolution is the epoch of exploring the structure and development of the near-surface part of the Earth's interior (to a depth of 10 km). This epoch, which started in the mid-19th century and continued almost until the end of the 20th century, can be called the epoch of geology, although a number of other sciences, namely, geophysics, geochemistry, mining, meteorology (physics of the atmosphere), and oceanology, formed in the 20th century. However, for a number of reasons, precisely the geology was a leading science of the Earth in the second epoch. [5] I am firmly convinced that the third epoch of the Earth science evolution started in the mid-1990s (tentatively, in 1995) and it should be called the epoch of geophysics and geoinformation science. [6] The main problem solved in this epoch is the exploring of the Earth as a whole: its structure down to the core, dynamic processes in the Earth's interiors, and its evolution during the last billion years. Moreover, in the third epoch the geoinformation science should be formed as an independent research branch, with its own subject of study and methods (this is discussed in detail below). [7] Two factors were of basic significance for the third epoch of Earth sciences: first, rapid development of computer technologies and, second, the creation of Earth science databases and information banks in the 1990s. [8] Because of volume limitations, this paper addresses only two topics: the essence of a new (developing) science, namely, geoinformatics, and the development of a basically new theory of interpretation of potential fields (gravity and magnetic anomalies). Section 1On the Essence of Geoinformatics[9] (1) As noted in Introduction, geoinformatics does not exist presently as an independent science because, although it has its own subject matter (The essence of this subject matter is systematization and accumulation of information that is related to the Earth and various geophysical fields and can be claimed in scientific research.), it does not have its own method as yet. As before, graphic (primarily cartographic) means of the geoinformation representation prevail. [10] (2) The aforesaid immediately raises the question of what is the essence of the method inherent in the geoinformation science. This question can be answered very briefly: the geoinformatics method consists in an analytical (naturally approximate) description of geoinformation. [11] However, such a brief answer will be unintelligible to the vast majority of readers. For this reason, the meaning of the answer is explained in detail below. [12] First of all, I present a rather evident classification of geoinformation. This classification includes the following four types of geoinformation: [13] (I) geoinformation defined on the set of land surface points; [14] (II) information related to the Earth's interiors; [15] (III) information related to water media (oceans and seas); and [16] (IV) information related to the Earth's atmosphere. [17] It is clear that the first two types of geoinformation are most significant. The following examples of these types can be presented. [18] Example 1. Information on the land topography of the Earth. [19] Example 2. Information on the Earth's gravity field (the complete field and its normal and anomalous components). [20] Example 3. Information on the magnetic field of the Earth (the complete field and its normal and anomalous components). [21] It is clear that these examples by no means exhaust the types of surface geoinformation. As regards examples of the second types of geoinformation, they are as follows. [22] Example 1. Information on the occurrence depths of interfaces between layers in the crust and primarily in the sedimentary cover of platform regions. [23] Example 2. Information on the spatial position of crustal faults. [24] (3) The essence of linear analytical approximations of surface geoinformation is explained fairly comprehensively in this and subsequent paragraphs of the section. [25] To begin with, note that geoinformation can be of three types: [26] (a) local (the sphericity of the Earth can be neglected); [27] (b) regional (geoinformation is specified on a fairly large area and the sphericity of the Earth should be taken into account); [28] (c) global (geoinformation is specified on the entire surface of the Earth). [29] Below I consider only the cases of surface geoinformation of the first type, implying that the Earth can be treated as a half-space bounded by a relatively small part of the Earth's surface. I address here the analytical description of the Earth's surface topography in a local variant. Let an orthogonal coordinate system
be introduced, with the 0 x3 axis directed upward (Figure 1) and let
be the elevations of the Earth's surface above the plane
which is regarded as the surface of a normal Earth. [30] It is clearly seen that the function h(x) can be treated as limiting values of a function harmonic in the half-space x3>0. It is in terms of this treatment that the sought linear analytical approximations can be constructed from given values:
where
and dhi are uncertainties in the given values hi. [31] There exist two main constructions, each involving a distribution of formal point sources on the plane
[32] First construction consists in an integral representation of the function h(x):
where
Below, this function is denoted as
[33] It is clear that the sought function in (1.7) is
r(x1, x2) (it should be found from the given values
hi,d,
1 [34] The function r(x1, x2) can be determined in terms of the constrained variational problem
Problem (1.10) is evidently solved by the method of Lagrangian multipliers, reducing it to a family of unconstrained variational problems of the type
where all values
li,
1 [35] Using the necessary (and, in this case, sufficient) criterion of an extremum and writing out the corresponding Euler equation [Buslaev, 1980; Kosha, 1979; Lavrentyev and Lyusternik, 1950], we obtain the following expression for the function r(x1, x2):
As can be shown (see Section 2 of this work), the vector
l with the components
li,
1
where the
N
and its elements are
It is easy to see that
The integrals on the right-hand sides of relations (1.15) have fairly complex analytical expressions and this is a significant drawback of the method described (which is a variant of the method of linear integral representations, described in a more general form in the next section of the paper). [36] We emphasize that the integrals determining the elements aij of the matrix A can be calculated by cubature formulas with a relatively small number of nodes (30 to 40); the integral
should be preliminarily transformed into the integral of the type
It is natural that we should actually find a stable approximate solution l of the SLAE
this solution should be consistent with a fairly large
volume of a priori information on the vectors
h (with the
components
hi,
1 [37] The main aspects of a new theory developed by the author for the determination of stable approximate solutions to SLAEs of form (1.19) are described in Section 2. [38] (4) Now, we characterize, in the local variant, the second (in a sense, more effective) method of constructing linear analytical approximations of the Earth's surface topography. [39] The function h(x) = h (x1, x2, 0) in this method is represented as
where we set
The points
coincide in coordinates (x(j)1, x(j)2) with the coordinates (x<undef>(j)1, x2<undef>(j)) of the points specifying the function h(x) = h(x1, x2, 0). In other words, a point source is located above each point (x1(j), x2(j), 0) at a depth H, and a set of such sources approximates the function h(x) = h(x1, x2, 0) treated in terms of limiting values of function harmonic in the exterior of the upper (relative to the normal surface of the Earth) half-space.
[40] It is clear that analytical expression (1.20) leads to a SLAE for
the determination of the coefficients
c(j),
1
hd is the
N -vector with the components
hi, d = hi + dhi,
c is the
N -vector of the sought
coefficients
c(j),
1
and consisting of the elements
Evidently, we have
and the matrix A strongly gravitates toward the banded type; i.e. the value | A | E is determined by the elements of the main diagonal and a small number of diagonals parallel to the main one. [41] Due to the above considerations, the determination of stable approximate solutions c to SLAEs of form (1.25) is not very difficult even in the case of large dimensions of the systems, namely,
[42] The following method is quite rational (effective). SLAE (1.23) is rewritten in the form
where A0 is a banded matrix with a band width 2m + 1,
DA is the matrix complementing the matrix A, and a>0 is a small number,
[43] The following iterative method is then used:
and
where we set
The quantities gn in (1.32) are found from the conditions
[44] Evidently, condition (1.32) can be written as
where q(n-1) and p(n-1) are SLAE solutions.
[45] The criterion for stopping this iterative process is defined by a priori information available for the vectors h and dh and is discussed in greater detail in the next section of the paper. [46] We believe that the readers are now convinced that the second method of constructing linear approximations of the Earth's surface topography (in a local variant) is much more effective than the first method for the following two reasons: [47] (a) the calculation of elements of the matrix A is much simpler; and [48] (b) the matrix A has specific properties facilitating the determination of the vector c. [49] (5) Now, after the problem of constructing linear analytical approximations of the Earth's surface topography (in a local variant) is considered in rather great detail, we address the problem of constructing linear analytical approximations of other types of geoinformation specified at the physical surface of the Earth (pertinent examples are given above in paragraph 2 of this section). Evidently, we should first consider this problem in a local variant. The only distinction of problems of constructing geoinformation on a physical surface S from the above problem of the topography approximation is that functions defined at the Earth's surface S should be treated as the limits of functions harmonic in the exterior of S (whereas the topography heights h(x) = h(x1, x2, 0) are approximated by functions harmonic above the plane x3 = 0). [50] Let
be the coordinates of points above the plane x3 = 0 at which approximate values of the function
are given; more specifically, these values can be written as
(We mean that u(x) is not a component of the gravitational or magnetic field but describes other geoinformation, for example, the intensity of the surface heat flow of the Earth.)
[51] Then, the coordinate system ( 0x1 x2 x3 ) is defined in
such a way that the axis
0x3 is directed downward and all points
x(i),
1
are the same as in points x(i), while the coordinates x(j)3 are equal in absolute value but opposite in sign to the corresponding coordinates of observation points (Figure 2):
Substituting the coordinates of observation points into
(1.41) and replacing
u(x(i)),
1
where
c is the
N -vector with the components
ci,
ud is the
N -vector with the components
ui,d,
and
A is the
N
and consisting of the elements
Evidently, we have
and hence
[52] It is clear that SLAE (1.44)-(1.46) is computationally rather simple, and the determination of its stable approximate solution c should not encounter great difficulties even in the case of SLAEs of large dimensions.
[54] Let we have n layers and, consequently, n+1 surfaces bounding these layers from above and from below. These boundaries can be described analytically if a sufficiently great number of points with known coordinates ( x1, x2, x3 ) in an a priori given coordinate system are present on each of the surfaces (Figure 3). [55] Note. If the boundaries of the layers lie at depths of about 1 km or more and the number of boreholes crossing these boundaries is fairly small (no more than 5-7), a sufficiently great number of interface points with known, albeit approximate, coordinates can be found solely geophysical, primarily seismic and electric, methods. This fact is a very weighty argument indicating that the new (third) epoch in the development of Earth sciences is an epoch of geophysics and geoinformatics. [56] An arbitrary boundary between layers G can be advantageously described in the following analytical form:
where Pm(x1, x2) and Qm(x1, x2) are algebraic polynomials of the given degree m,
All coefficients kp, q and cp, q in (1.50) and (1.51) should be found from the given values
The determination of the coefficients kp, q and cp, q evidently reduces to the solution of the SLAE derived from the relations
It is clear that the degrees m of the polynomial Pm(x1, x2) and Qm(x1, x2) should be small. The most important is the case of
In this case, the number of coefficients to be determined is 19; therefore, it is convenient to have the same (even greater) number of the points x(i) = (x(i)1, x2(i), x3(i)) with known coordinates on a boundary G between layers. It is clear that this condition is valid in the majority of cases. [57] (7) Finally, we have to consider (as briefly as possible) the case when the surfaces to be analytically approximated are fault planes and boundaries of geological bodies. [58] It is appropriate to specify analytical approximations of G in the following general form:
where F(x1, x2, x3) is either an algebraic polynomial or any other function including linear coefficients to be determined from coordinates of the points given on the surface G. It is clear that these coefficients should be found from a system of linear equations of the form
The system can be either normally determined (if the number of the known points on the surface G is the same as the number the coefficients to be found) or overdetermined (if the former is greater than the latter). The variant of an overdetermined SLAE is evidently preferable. However, the required number of points on G can actually be obtained from detailed geophysical investigations, preferably by a complex of methods. [59] This is an additional argument indicating that the new (third) epoch in the development of Earth sciences is an epoch of geophysics and geoinformatics. Section 2Development of a Basically New Theory of Interpretation of Potential Fields (Gravity and Magnetic Anomalies)[60] (1) The study of anomalous gravitational and magnetic fields began in the 19th century in Russia (examples are the studies of the Moscow gravity anomaly and the Kursk magnetic anomaly [Gamburtsev, 1960a; Ilyina, 1983]. However, gravimetry and magnetometry on the whole and methods of potential field interpretation in particular were most intensely developed in the 20th century. [61] Following the well-known concept of Thomas Kun [6], two paradigms and five periods are recognizable in the development of the theory of potential field interpretation in the 20th century. These paradigms, characterized below, are (1) paradigm of manual computations and (2) paradigm of the early computer epoch. [62] As regards periods in the development of the theory of potential field interpretation in the 20th century, they are as follows: [63] I period (1900-1919), pre-paradigm period; [64] II period (1920-1939), buildup period of the epoch of manual computations; [65] III period (1940-1959), period of "normal science" within the framework of the manual computation paradigm; [66] IV period (1960-1985), buildup period of the paradigm of the early computer epoch (It is clear that the boundaries of the periods are tentative but cannot be shifted by more than two to three years forward or backward.) [67] V period (1986-2000), period of "normal science" within the framework of the paradigm of the early computer epoch. It is the author's strong conviction that third paradigm, that of the mature computer epoch, has started to develop since 2001. [68] (2) The main aspects of the basically new theory of interpretation of potential fields (gravity and magnetic anomalies) are as follows. [69] I. Leading role of the general methodology of interpretation. Within the framework of this general methodology, the concept of method-generating ideas is of the greatest significance. These are the ideas of (a) approximation, (b) linearization, and (c) optimization. [70] In particular, it is emphasized in general methodology that the main equations describing field potentials and their first derivatives with respect to Cartesian coordinates are linear differential equations of Laplace and Poisson (in the exterior of the carrier of field sources and inside the carrier, respectively). [71] II. Importance of the "geophysical dialect" of the language of mathematics. The essence of the "geophysical dialect" concept is that mathematical constructions accounting for specific features and problem formulations should be used in all geophysical methods (first of all, in gravimetry and magnetometry). [72] III. Methods of interpretation of potential fields should be based on two basic mathematical apparatuses: [73] (a) classical continual theories of potential fields, and [74] (b) discrete theories of potential fields proposed by the author. [75] Together with his disciples, the author has published a large series of studies on the application of methods of discrete theories of gravitational and magnetic fields to the solution of interpretation problems. [76] IV. Reduction of any metrological and interpretation problems of gravimetry and magnetometry to the determination of stable approximate solutions of SLAEs of the form
the vectors on the right-hand side are approximately specified data of observations of anomalous potential fields under study. [77] We should emphasize that, within the framework of the new theory of interpretation of potential fields (gravity and magnetic anomalies), the solution of inverse problems is also reduced to the solution (generally iterative) of SLAEs. As shown in detail below, the finite element description of the geological medium under study (within the framework of classical continual theories) is of paramount importance in the new theory. [78] V. Use of a basically new information base in gravimetry and magnetometry. The essence of the new information base is the use of linear analytical approximations of observed anomalous fields. More specifically, we mean the replacement of maps (albeit digital) of observed fields by their linear analytical approximations, which are preferable in a number of aspects. [79] Three main representations used for constructing linear analytical approximations of potential fields in local and regional variants are as follows: [80] (1) set of point source fields; [81] (2) segments of series in spherical functions; and [82] (3) linear integral representations whose integrands are a linear combination of given functions. [83] VI. Use of a symmetrized (relative to the vertical axis 0z ) field. This approach enhances the efficiency of methods of theories of gravitational and magnetic fields. [84] VII. Basically new theory of the determination of stable approximate solutions of SLAEs (2.1) that is fully adequate to the conditions and demands of geophysical (first of all, gravimetric and magnetometric) practice. This theory [Strakhov, 1997a, 1997b, 1997c; and others] has been developed because the classical theory of regularization of SLAEs (2.1) is inadequate to the conditions and demands of geophysical practice. (We mean the theory of the SLAE regularization, which has been created in fundamental works by A. N. Tikhonov, M. M. Lavrentyev, and V. K. Ivanov, as well as by their numerous disciples and followers). [85] (3) Now, after formulating the main concepts of the theory of potential fields (gravity and magnetic anomalies) in terms of the third paradigm, it is appropriate to demonstrate its basic distinction from theories that have been developed within the framework of the first and second paradigms. [86] In accordance with the name itself of the first paradigm (the epoch of manual computations), its theory of interpretation of gravity and magnetic anomalies was primitive. This is evident from the following considerations. [87] (A) The theory was oriented toward plane (2-D) problems. [88] (B) Various manual tools such as master plots, nomograms, and graphic means played the decisive role. [89] (C) Forward problems were solved only for very simple bodies having a uniform density or magnetization. [90] (D) The reduction of interpretation problems to the SLAE solution was not used at all. [91] Now, we characterize the essence of the second paradigm in the theory of potential field interpretation. Here we should note first of all the name of the paradigm (the early computer epoch) and the fairly long period of its development, which is related to gradual (from the modern standpoint, very slow) progress in computer technologies. [92] We should emphasize already here that only one basic mathematical apparatus used in the second paradigm was the apparatus of classical continual theories of gravitational and magnetic fields. The following characteristic features are inherent in the potential field interpretation theory of the second paradigm. [93] (A) In addition to methods of anomalous field interpretation in terms of 2-D (plane) problems, interpretation methods started to be developed in terms of 3-D problems, mainly for the cases when the sphericity of the Earth can be neglected. [94] (B) The fitting method was generally used for solving inverse problems. In this respect, methods of solving forward problems were intensely developed in both 2-D and 3-D variants. [95] (C) Much progress was made in studying problems of equivalence and uniqueness in gravity and magnetic data inversion in both 2-D and 3-D variants. [96] (D) Methods of anomalous potential field interpretation based on the apparatus of spectral representations of potential field elements were largely developed [Gamburtsev, 1960b; Gladkiy and Serkerov, 1974; Serkerov, 1991], (V. N. Strakhov and A. I. Luchitskiy, preprint, 1980). [97] (E) Along with methods of solving forward and inverse problems, much attention was given to potential field transformation methods, primarily, the methods of analytical continuation of potential fields into the lower half-space and the related determination of singular fields of the fields [Strakhov, 1984]. [98] (F) The theory of 2-D potential field interpretation was basically transformed: its elaboration was now based on the theory of functions of complex variables [Shalaev, 1960; Tsirulskiy, 1990]. [99] (G) The theory of interpretation of potential fields started to use methods of the theory of ill-posed problems, developed by A. N. Tikhonov, M. M. Lavrentyev, and V. K. Ivanov, as well as by their numerous disciples and followers. [100] On the whole, we can state that the potential field interpretation theory of 1985 basically differed from that existing in 1959, but its development in Russia in the fifth period (1986-2000) dramatically slowed down, first of all, because of the economic depression (However, there is no doubt that a stereotype of thinking formed (tentatively by 1990) in the interpretation theory of potential fields played also a certain role in this stagnation tendency.). [101] (4) Now, we address the essence of the basically new theory of potential field interpretation developed within the framework of the third paradigm. Below, we describe, first, the method of linear integral representations proposed by the author in the mid-1990s and, second, the method of gravity data inversion based on finite element descriptions of the geological medium studied. For simplicity (and this is the only reason), our description is given with reference to 2-D (plane) problems. [102] (5) Now, we address the essence of the linear integral representation method, generalizing the classical method of linear integral equations of first kind. [103] Let the N numbers
be given; here, fi are useful signals and dfi are their determination uncertainties. We assume that all functions fi can be expressed analytically as
where
Mr are given connected point sets in
R2 (2-D
Euclidean space);
mr(x) are measures on the sets
Mr;
Qr(i)(x),
1 [104] This problem is solved on the basis of the formulation and solution of the following constrained variational problem:
where all functions
pr2(x) are given and account for
a priori information on the sought functions
rr(x). In
particular, if the sets
Mr (all or some of them) have boundaries
should be valid, the functions pr2(x) must be such that
[105] It is clear that constrained variational problem (2.4) can be solved most simply by the classical method of Lagrangian multipliers, i.e. in terms of the family of unconstrained variational problems
so that the numbers li in the functional to be minimized are the Lagrangian multipliers accounting for the conditions inherent in constrained variational problem (2.4).
[106] Using the classical necessary (and in this case sufficient)
minimization conditions, i.e. constructing Euler equations
[Buslaev, 1980;
Gamburtsev, 1960a, 1960b],
we obtain explicit analytical expression of the sought functions
rr(x),
1
The numbers
li,
1
where
f is the
N -vector with the components
fi,
1
and its elements are
However, since all components
fi,
1
Thus, the problem of determining
R functions
rr(x),
1 [107] (6) There exists another, basically important variant of the linear integral representation method. This variant reduces to the fact that the functions
are a priori known and the main problem is their refinement. In this variant, constrained variational problem (2.4) is replaced by the constrained variational problem
where
and
[108] In this case as well, the method of Lagrangian multipliers is used to reduce problem (2.14) to the family of unconstrained variational problems
This readily yields approximate (but sufficiently accurate) expressions for the sought functions Drr(x):
Based on relations (2.17) and conditions in (2.13), we
find that the
N -vector
g with the components
gi,
1
where the
N
and its elements are
However, because the vector Df is unknown, and the known vector is
we should actually find a stable approximate solution of the SLAE
Finally, we should only mention the validity of the equality
where A is the matrix in SLAEs (2.9) and (2.12). [109] (7) Here, we present a very brief description of one more variant of the linear integral representation method that is important for joint interpretation of gravity and seismic data.
[110] This variant also implies that zeroth approximations
rr(0)(x),
1
The values cr are known,
It is clear that problem (2.25) should also be solved by the method of Lagrangian multipliers. [111] Omitting details sufficiently clear for the reader, the solution of the constrained variational problem is obtained in the form
where
The values
hr,
1
The technique of reducing the problem of determination of N + R values
to the determination of a stable approximate solution of a certain SLAE should quite clear to the reader, so that we do not write out here the corresponding SLAE.
[112] In conclusion, we note that there exists a number of alternative
formulations of the problem of determining the functions
rr(x),
1 [113] Below, we present three examples demonstrating the application of the linear integral representation method to the solution of gravity problems in the 2-D (plane) variant because 2-D problems admit clear and simple 2-D graphics. [114] (8) First example illustrates the problem of analytical continuation of 2-D potential fields. Let at points
of the function Dg(x, z)
where Va(x, z) is the potential of the anomalous gravitational field produced by field sources in the lower half-plane (Figure 4). We emphasize that the 0z axis is directed downward. Thus, we have
and
whereas
dDgi,
1 [115] Evidently, it is possible (and appropriate) to represent the analytical approximation of the function Dg(x, z) as the following Fourier integral:
where w max is a fairly large number; to choose the latter, a priori information on the depth to the upper edge of anomalous masses and the values
are used. [116] The function Dg(x, z) is analytically represented in the standard form used in the method of linear integral representations. We have
and
[117] Evidently, we have in this case:
[118] The SLAE used for the determination of the vector
l (with the components
li,
1
where the
N -vector
Dgd has the components
Dgi, d,
1
and its elements are expressed analytically as
It is known that
and, therefore, we obtain the explicit analytical expressions for the elements aij of the matrix A
[119] Thus, all elements aij of the matrix A have rather simple analytical expressions. A reasonable value for w max is
where
[120] A stable approximate solution l of SLAE (2.42)-(2.43) can be found by well-known methods.
[121] After the vector
l is found, analytical expressions of
the functions
A(w) and
B(w) are also determined from
(2.42), where
li should be replaced by
li,
1
can also be found with the required accuracy at any level
in the upper (z < 0) and upper (z > 0) half-planes. As is well known, this also implies the possibility of finding singular points of the analytical continuation of field elements into the lower half-space, thereby gaining information carriers of anomalous masses that generated the anomalous gravitational field.
[123] Evidently, the function Dg(x, z) can be represented by the analytical expression
where rk(x, z) the densities of anomalous masses multiplied by a known constant. Further, the following conditions should be valid at the boundaries D Dk of the regions Dk:
Let the following values be known:
where
Here S is the earth-air interface (see Figure 5). [124] The method of linear integral representations in the given problem consists in the solution of the constrained variational problem
The functions p2k(x, z) in the minimization functional satisfy the conditions
[125] Note. Evidently, the functions p2k(x, z) can be taken in the form
where jk(x, z) are functions conformally mapping the regions Dk onto the unit circle. [126] The remaining calculations should be clear to the reader. Namely, we obtain the following expressions for the functions rk(x,z):
Now, we have to find a stable approximate solution of the SLAE
where the matrix A possesses the property
and consists of the elements
[127] It is clear that the integrals in (2.62) should be calculated by the
corresponding cubature formulas, which is a rather complex
procedure. However, "the game is worth the candle" because, if the
vector
l (a stable approximate solution of SLAE (2.60))
is found, all functions
rk(x, z),
(x, z)
[128] Further, once the functions
rk(x, z) are determined, the
fields
Dgk(x, z) of anomalous masses in the regions
Dk can
also be found (approximately but with the required accuracy), which
solves the problem of field separation. Moreover, the found
functions
rk(x, z),
(x, z) [129] Note. Evidently, the simplest cases are those in which the regions Dk are circles or ellipses because the functions jk(x, z) are then expressed by simple formulas.
![]() [131] Thus, let approximate values of the function Dg(x, z) given at points (x(i), z(i)) of the earth-air interface S be
these values are due to density distributions
rr(x,z) in the a priori given regions (blocks)
Dr,
r = 1, 2, [132] A constrained extremal problem is formulated as
where all functions vr(x, z) are known and
[133] Problem (2.64)-(2.65) is the constrained variational problem considered in paragraph 5, where it is shown to reduce to the solution of a certain SLAE.
In solving the gravity inversion problem, a constant density equal or unequal to zero is ascribed to each of the squares, and these density values are to be determined from an external anomalous field specified as a set of the values
We assume that the number N is sufficiently great, namely,
[135] The finite element description of the geological medium should meet the following conditions: [136] (a) each carrier of sources should be covered by a sufficiently large number of squares (finite elements) of the same side h;
[137] (b) the total number of squares covering all carriers of field
sources should not overly large, namely, it should be two to three
times smaller than the total number of the field values
Dg(x, z) calculated at all horizons
Gj,
1
[138] Note. In order to determine the field values at points of the
lines
Gj,
1 [139] We should specially emphasize that, if a square with a small side h is centered at a point (x, z), then even at the distance l
from this point, the approximate formula
ensures the required high accuracy (here, r is the density of masses in the square and G is the gravitational constant). Therefore, enumerating squares containing anomalous masses through the index k,
and denoting their masses as
and their centers as
the field Dg(x, z) of this set of squares can be described by the approximate expression
which ensures the required high accuracy. This immediately
implies that the masses
mk,
1
where
m is the
N1 -vector with the components
mk,
k=1, 2,
and
A is the
N
where
are coordinates of points (of the lines
Gj,
1 [140] Evidently, it makes no sense (for many reasons) to find exact solution of SLAE (2.77). Its approximate solution should be found under the condition
where D2 is an a priori given quantity. For this purpose, it is appropriate to use iterative methods. [141] (12) In this paragraph, we consider (perforce, very briefly) the main concepts of the new theory of stable approximate solutions x to SLAEs of the form
where
A is an
N [142] The vector x is a stable approximate solution of SLAE (2.82) if the following approximate relations are satisfied with sufficiently high accuracy:
[143] Evidently, the vector x satisfying relations (2.83) can be found only if a sufficiently large volume of a priori information on the vectors df and f is available. [144] First type of the a priori information is the knowledge of the constants d min2 and d2 max in the inequalities
Moreover, the ratio
suggests the appropriateness of considering two situations. [145] First situation, which simpler and more widespread, is specified as
[146] Second situation, more complex and less widespread, is characterized by the interval
[147] Second type of the a priori information is the knowledge that the vectors of the useful signal f and noise df are mutually orthogonal:
[148] Third type of the a priori information (being of fundamental importance!) is the knowledge that the noise vector df is random and homogeneous. The latter implies that
[149] (a) the numbers of positive and negative components
dfi,
1
are nearly the same;
[150] (b) all components of
dfi,
1
[151] (c) all components
dfi,
1 [152] Fourth type of the a priori information is the knowledge of the functional
establishing the preference relation on the set of
approximate solutions of SLAE (2.82); here,
G is an
N [153] Comment. The functional W(x) establishes a preference relation on the set of approximate solutions of SLAE (2.82) if, for any two approximate solutions x(1) and x(2) satisfying the condition
the inequality
implies that the approximate solution x(2) is preferable.
[154] In the first situation ( 1 < g2
[155] Thus, we have two different approaches to the determination of the vector x. [156] In the first approach (the only approach used at present), the vector x is first found, after which the vectors u and z are determined:
[157] In the second approach, the vector
is first determined, and the vector x is then found as a sufficiently accurate approximate solution of the SLAE
As shown below, the second approach is preferable in
the case of SLAE (2.82) of a large dimension (tentatively,
min(N, M) [158] The point is that the vector u can (and must!) be found from the solution of the constrained extremal problem
Furthermore, it is evident (and basically important!) that the conditions
can be replaced by the conditions
Therefore, constrained extremal problem (2.97) is replaced by the following one:
Using the method of Lagrangian multipliers, constrained extremal problem (2.97) corresponds to the family of unconstrained extremal problems
where
are the Lagrangian multipliers accounting for the equality conditions to problem (2.100). [159] The solution ul, m of problem (2.101) evidently satisfies the SLAE
[160] Since
GG T is an
N
where
and
Here
EN is the unit
N [161] The value l is always sufficiently large,
Consequently, the matrix Sl is sufficiently well-conditioned at all l and the solution of constrained extremal problem (2.100) should not involve great difficulties. [162] It is clear that, using the classical compact Gaussian scheme, the matrix Sl can be represented as the product of triangular matrices
where Ll is the lower triangular matrix, with all diagonal elements being equal to unity, and Rl is the upper triangular matrix. [163] Thus, the solution of SLAE (2.104) reduces to the successive solution of two SLAEs with the triangular matrices:
[164] Therefore, it is sufficient to find a pair (l, m) at which the conditions (equalities) to problem (2.102) are satisfied (approximately but with required accuracy). [165] At first glance, this problem appears to very cumbersome. However, taking into account the structure of the vector on the right-hand side of SLAE (2.104), it can be significantly simplified. Namely, the vector ul, m can be represented in the form
where ul, 0 is the SLAE solution at m = 0 and ul, 1 is the SLAE solution in the case if the vector on the right-hand side has the form
where
Now, the value ml in (2.112) is defined in such a way that it ensures the validity of the second condition in constrained extremal problem (2.104) at any l. As is easily seen, this value is
[166] Thus, only the value l is to be varied (for finding a solution that makes both conditions in (2.100) valid). Evidently, if computations are made on PC clusters or multiprocessor systems, parallel computations at different values of l,
are possible, which saves the computing time very significantly. [167] The values of l can be divided into the following two groups: [168] (a) a localizing group that defines the required range of l values and [169] (b) a solving group that is generated from the results of computations of the first group and ensures the determination of the sought (final) value. [170] Now, we discuss the technique of determining an approximate (but sufficiently accurate) solution of the SLAE
[171] It is clear that the most rational approach is to use new iterative methods. The criterion of terminating the process of successive iterations
should be taken in the form
where Crit is an a priori given number related to the value
Namely, it is appropriate to take
[172] (13) Concluding Section 2 of this work, we describe a basically new method of finding the vector x within the framework of the first approach according to which the problem of filtering of the given vector fd is not solved and the vector x is found directly from the vector fd with regard for the available a priori information. [173] The new method is based on the formulation of the following constrained extremal problem:
where the quantity p2 is given and
[174] We emphasize that the formulation of this constrained extremal problem is novel with respect to [175] (a) the minimization functional, and [176] (b) the equality conditions. [177] Note. The conditions in problem (2.121) are equivalent to the conditions
but are much simpler analytically. [178] Constrained extremal problem (2.121) should again be solved by the method of Lagrangian multipliers, i.e. by considering the family of the unconstrained extremal problems
where
are the Lagrangian multipliers accounting for the equality conditions in problem (2.121). [179] The solution of problem (2.124) satisfies the SLAE
where
and, because
G is an
N
where
and
Here,
EN+m and
EM are unit matrices of the sizes
(N+m) [180] As in the case of SLAE (2.128), the vector xl, m can be represented in the form
where xl, 1 and xl, 0 are determined through the respective SLAE solutions
It is clear that, for any given l, the matrix Wl is decomposed (according to the classical Gaussian scheme) into the product of triangular matrices
where
Ll is the lower triangular matrix of the size
(N + m + M) [181] In conclusion, we emphasize that, if the values of l are not overly large, matrix Wl (2.131) is sufficiently well-conditioned (although to a lesser degree as compared with matrix Sl (2.107) used in the filtering problem). 3. Conclusion[182] To sum up, the following seven results of this work are most significant. [183] (1) Undoubtedly, the most important result is the substantiation of the concept according to which three epochs are distinguished in the evolution of Earth sciences and the third epoch of geophysics and geoinformatics started to develop in 1995. [184] (2) The paper described the subject and method of the new science of geoinformatics. [185] (3) A new method is described for the representation of geoinformation given at the Earth's surface; in terms of this method, the information is represented through linear analytical approximations, based on the treatment of functions defined at the Earth's surface as the limiting values of functions harmonic in the exterior of the Earth. [186] (4) A new method is proposed for the description of the Earth's interiors in the form of linear analytical approximations of surfaces of geological bodies and faults. [187] (5) The paper presented the main concepts of a basically new theory of interpretation of potential fields (gravity and magnetic anomalies) that defines the third paradigm of the potential field interpretation. [188] (6) Two main techniques of reducing gravity problems to the solution of a system of linear analytical equations (SLAE) on the basis of [189] (a) the method of linear integral representations (in several variants) proposed by the author, and [190] (b) a finite element description of a volume of the geological medium under study. [191] We should emphasize that gravity inversion problems are also reduced the SLAE solution. [192] (7) The very important result of this paper is a fairly detailed description of the basically new theory of finding approximate solutions x to the SLAE
and two new methods of determination of the vector x with reference to the case
where d2 min and d2 max are the constants in the inequalities
[193] Both methods are based on [194] (i) the formulation of new constrained extremal problems, and [195] (ii) the reduction of these problems to SLAEs of a new type, namely, to extended SLAEs of second kind. [196] However, it is clear that a great deal of effort must be made to develop computer technologies implementing all theoretical ideas proposed in the paper and to test these technologies in a large number of model experiments and on real data. ReferencesBuslaev, V. S. (1980), Variational Calculus (in Russian), 285 pp., LGU, Leningrad. Gamburtsev, G. A. (1960a), On the origin of the Kursk gravity and magnetic anomaly: Attraction by ridges (in Russian), 72 pp., Akademizdat, Moscow. Gamburtsev, G. A. (1960b), Geomagnetic interpretation of gravity and magnetic observations with the use of mechanical computing devices, 1, Selected Papers (in Russian), 81 pp., Akademizdat, Moscow. Gladkiy, K. V., and S. A. Serkerov (1974), Fourier Transformations and Their Applications in Gravity and Magnetic Surveys (in Russian), 72 pp., MINXiGT, Moscow. Ilyina, T. D. (1983), Formation of the Soviet School of Geophysical Exploration (in Russian), 215 pp., Nauka, Moscow. Kosha, A. (1983), Variational Calculus (in Russian), 279 pp., Vysshaya Shkola, Moscow. Lavrentyev, M. A., and L. A. Lyusternik (1950), Course of Variational Calculus (in Russian), 296 pp., Gostoptekhizdat, Moscow. Serkerov, S. A. (1991), Spectral Analysis in Gravity and Magnetic Surveys (in Russian), 279 pp., Nedra, Moscow. Shalaev, S. V. (1960), Application of functions of complex variables to geological interpretation of gravity and magnetic anomalies, Vopr. Razved. Geofiz. (in Russian), 1, 3. Strakhov, V. N. (1984), Methods of Interpretation of Gravity and Magnetic Anomalies: A Handbook of the Special Course (in Russian), 72 pp., PGU, Perm. Strakhov, V. N. (1997a), General theory of stable approximate solutions to systems of linear algebraic equations with approximately given right-hand sides and matrices related to geophysical problems, in: Theoretical and Applied Problems of Geological Interpretation of Gravitational, Magnetic, and Electric Fields (in Russian), p. 38, OIFZ RAN, Moscow. Strakhov, V. N. (1997b), Mathematical apparatus used in the construction of algorithms of finding stable approximate solutions to systems of linear algebraic equations in gravity and magnetic problems, in: Theoretical and Applied Problems of Geological Interpretation of Gravitational, Magnetic, and Electric Fields (in Russian), p. 43, OIFZ RAN, Moscow. Strakhov, V. N. (1997c), Extremal problems, nonparametric regularization, and filtering in the theory of stable approximate solutions to systems of linear algebraic equations with approximately given right-hand sides and matrices, in: Theoretical and Applied Problems of Geological Interpretation of Gravitational, Magnetic, and Electric Fields (in Russian), p. 76, OIFZ RAN, Moscow. Tsirulskiy, A. V. (1990), Functions of Complex Variables in the Theory and Methods of Interpretation of Geophysical Potential Fields (in Russian), 132 pp., URO AN SSSR, Sverdlovsk. Received 21 February 2007; revised 30 February 2007; accepted 22 March 2007; published 10 August 2007. Keywords: epochs in earth sciences, geoinformatics, gravity field analysis. Index Terms: 0520 Computational Geophysics: Data analysis: algorithms and implementation; 0903 Exploration Geophysics: Computational methods: potential fields; 0920 Exploration Geophysics: Gravity methods; 1217 Geodesy and Gravity: Time variable gravity; 3200 Mathematical Geophysics. ![]() Citation: 2007), Change of epochs in Earth sciences, Russ. J. Earth Sci., 9, ES1001, doi:10.2205/2007ES000217. (Copyright 2007 by the Russian Journal of Earth SciencesPowered by TeXWeb (Win32, v.2.0). |