RUSSIAN JOURNAL OF EARTH SCIENCES VOL. 9, ES1001, doi:10.2205/2007ES000217, 2007

Section 2

Development 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

eq057.gif(2.1)

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

eq058.gif(2.2)

be given; here, fi are useful signals and dfi are their determination uncertainties. We assume that all functions fi can be expressed analytically as

eq059.gif(2.3)

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 le r le R, 1 le i le N, are given functions defined on the corresponding sets Mr; and rr(x), 1 le r le R, are functions to be determined (approximately but with sufficient accuracy) from the given values fi, d, 1le ile N.

[104]  This problem is solved on the basis of the formulation and solution of the following constrained variational problem:

eq060.gif(2.4)

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 partial Mr at which the conditions

eq061.gif(2.5)

should be valid, the functions pr2(x) must be such that

eq062.gif(2.6)

[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

eq063.gif(2.7)

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 le r le R, through the Lagrangian multipliers:

eq064.gif(2.8)

The numbers li, 1 le i le N are components of an N -vector called the Lagrange vector and, using conditions in (2.4), we obtain the following SLAE for the determination of these components:

eq065.gif(2.9)

where f is the N -vector with the components fi, 1 le i le N. As regards the Ntimes N matrix A, it possesses the property

eq066.gif(2.10)

and its elements are

eq067.gif(2.11)

However, since all components fi, 1 le i le N, of the vector f are unknown and only the components fi, d of the N -vector fd are known, we should actually find a stable approximate solution l of the SLAE

eq068.gif(2.12)

Thus, the problem of determining R functions rr(x), 1 le r le R, is reduced to the determination of a stable approximate solution of SLAE (2.12) alone.

[107]  (6) There exists another, basically important variant of the linear integral representation method. This variant reduces to the fact that the functions

eq069.gif(2.13)

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

eq070.gif(2.14)

where

eq071.gif(2.15)

and

eq072.gif(2.16)

[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

eq073.gif(2.17)

This readily yields approximate (but sufficiently accurate) expressions for the sought functions Drr(x):

eq074.gif(2.18)

Based on relations (2.17) and conditions in (2.13), we find that the N -vector g with the components gi, 1 le i le N, satisfies the SLAE

eq075.gif(2.19)

where the Ntimes N matrix B possesses the property

eq076.gif(2.20)

and its elements are

eq077.gif(2.21)

However, because the vector Df is unknown, and the known vector is

eq078.gif(2.22)

we should actually find a stable approximate solution of the SLAE

eq079.gif(2.23)

Finally, we should only mention the validity of the equality

eq080.gif(2.24)

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 le r le R, of the functions rr(x) are known, but the constrained extremal problem is formulated with the use of the measure of correlative closeness of the sought functions rr(x) to the known functions rr(0)(x), 1 le r le R. Namely, the problem of the determination of the functions rr(x) is formulated as

eq081.gif(2.25)

The values cr are known,

eq082.gif(2.26)

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

eq083.gif(2.27)

where

eq084.gif(2.28)

The values hr, 1 le r le R, like the values si, 1 le i le N, are unknown and should be determined from the values fi, d, 1 le i le N. It is evident that, for any gravity problems, we have

eq085.gif(2.29)

The technique of reducing the problem of determination of N + R values

eq086.gif(2.30)

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 le r le R, but they cannot be described here.

[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

eq087.gif(2.31)

2007ES000217-fig04
Figure 4
of the earth-air interface (Figure 4) be known the values

eq088.gif(2.32)

of the function Dg(x, z)

eq089.gif(2.33)

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

eq090.gif(2.34)

and

eq091.gif(2.35)

whereas dDgi, 1 le i le N, are as before the determination uncertainties of the values Dgi.

[115]  Evidently, it is possible (and appropriate) to represent the analytical approximation of the function Dg(x, z) as the following Fourier integral:

eq092.gif(2.36)

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

eq093.gif(2.37)

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

eq094.gif(2.38)

eq095.gif(2.39)

eq096.gif(2.40)

and

eq097.gif(2.41)

[117]  Evidently, we have in this case:

eq098.gif(2.42)

[118]  The SLAE used for the determination of the vector l (with the components li, 1 le i le N ) is

eq099.gif(2.43)

where the N -vector Dgd has the components Dgi, d, 1 le i le N, the Ntimes N matrix A possesses the property

eq100.gif(2.44)

and its elements are expressed analytically as

eq101.gif(2.45)

It is known that

eq102.gif(2.46)

and, therefore, we obtain the explicit analytical expressions for the elements aij of the matrix A

eq103.gif(2.47)

[119]  Thus, all elements aij of the matrix A have rather simple analytical expressions. A reasonable value for w max is

eq104.gif(2.48)

where

eq105.gif(2.49)

[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 le i le N. However, this also implies that the functions Dg(x, z) and

eq106.gif(2.50)

can also be found with the required accuracy at any level

eq107.gif(2.51)

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.

2007ES000217-fig05
Figure 5
[122]  (9) Second example relates to problems of separation of 2-D anomalous gravitational fields and determination of integral characteristics of partial carriers of anomalous masses. Let the field Dg(x) be generated by a certain number of local carriers of anomalous masses, for example, by three carriers Dk, k = 1, 2, 3 (Figure 5). We assume now that finite regions Dk, k = 1, 2, 3, undoubtedly containing the anomalous mass carriers Dk (Figure 5) are known.

[123]  Evidently, the function Dg(x, z) can be represented by the analytical expression

eq108.gif(2.52)

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:

eq109.gif(2.53)

Let the following values be known:

eq110.gif(2.54)

where

eq111.gif(2.55)

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

eq112.gif(2.56)

The functions p2k(x, z) in the minimization functional satisfy the conditions

eq113.gif(2.57)

[125]  Note. Evidently, the functions p2k(x, z) can be taken in the form

eq114.gif(2.58)

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):

eq115.gif(2.59)

Now, we have to find a stable approximate solution of the SLAE

eq116.gif(2.60)

where the matrix A possesses the property

eq117.gif(2.61)

and consists of the elements

eq118.gif(2.62)

[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) in Dk, can also be determined from (2.59).

[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)in Dk, enable rather simple determination of harmonic moments of masses with respect to certain given points inside the regions Dk.

[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.

2007ES000217-fig06
Figure 6
[130]  (10) Third example is undoubtedly the most important. In this case, the density distribution is to be found in a subhorizontal layer Dk bounded by known curves K1 and K2 (Figure 6), provided that the P wave velocity distribution in Dk is known and it is sufficiently strongly correlated with the density distribution. Moreover, it is a priori known that the degree of correlation between the P wave velocity v(x, z) and density r(x, z) differs, to an extent, in different blocks Dr, r=1, 2,ldots, R.

[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

eq119.gif(2.63)

these values are due to density distributions rr(x,z) in the a priori given regions (blocks) Dr, r = 1, 2,ldots, R.

[132]  A constrained extremal problem is formulated as

eq120.gif(2.64)

where all functions vr(x, z) are known and

eq121.gif(2.65)

[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.

2007ES000217-fig07
Figure 7
[134]  (11) Now we address the new method of gravity data inversion based on a finite element description of the studied part of the geological medium. As before, we consider the 2-D (plane) variant of the problem (see Figure 7). The finite element description of the studied (2-D) portion of the geological medium consists in its representation as a set of small squares with the side h. Depending on the concrete problem, the latter varies within the range

eq122.gif(2.66)

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

eq123.gif(2.67)

We assume that the number N is sufficiently great, namely,

eq124.gif(2.68)

2007ES000217-fig08
Figure 8
Evidently, formulas (2.65)-(2.68) can yield field values Dg(x, z) on a number of horizontal lines lying well above the sought carriers of anomalous masses (see Figure 8). The number of points spaced at a step Dx should be fairly large at each level z = const: more specifically, it should exceed the number of initial points at the earth-air interface S. The number n of horizontal lines Gj should lie within the (somewhat arbitrary) range

eq125.gif(2.69)

[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 le j le n (from the values Dgi, d = Dgi + dDgi at the earth-air interface points (x(i), z(i)) in S ).

[138]  Note. In order to determine the field values at points of the lines Gj, 1 le j le n, from the field values Dgi, d given at points on S, a linear analytical approximation of the field Dg(x, z) should be found in the form of a Fourier integral with finite limits, after which the values of the integral at the given points can be determined (it can be shown that the values of Fourier integrals are expressed in elementary functions).

[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

eq126.gif(2.70)

from this point, the approximate formula

eq127.gif(2.71)

eq128.gif(2.72)

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,

eq129.gif(2.73)

and denoting their masses as

eq130.gif(2.74)

and their centers as

eq131.gif(2.75)

the field Dg(x, z) of this set of squares can be described by the approximate expression

eq132.gif(2.76)

which ensures the required high accuracy. This immediately implies that the masses mk, 1 le k le N1 should be found from the solution of the SLAE

eq133.gif(2.77)

where m is the N1 -vector with the components mk, k=1, 2,ldots, N1; Dg is the N -vector with the components

eq134.gif(2.78)

and A is the Ntimes N1 matrix with the elements

eq135.gif(2.79)

where

eq136.gif(2.80)

are coordinates of points (of the lines Gj, 1 le j le n ) at which values of the function Dg(x, z) are known.

[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

eq137.gif(2.81)

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

eq138.gif(2.82)

where A is an Ntimes M matrix (in general, N M ), with N and M being the numbers of its rows and columns, respectively (it is clear that N is the number of the SLAE equations and M is the number of components of the sought vector x ).

[142]  The vector x is a stable approximate solution of SLAE (2.82) if the following approximate relations are satisfied with sufficiently high accuracy:

eq139.gif(2.83)

[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

eq140.gif(2.84)

Moreover, the ratio

eq141.gif(2.85)

suggests the appropriateness of considering two situations.

[145]  First situation, which simpler and more widespread, is specified as

eq142.gif(2.86)

[146]  Second situation, more complex and less widespread, is characterized by the interval

eq143.gif(2.87)

[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:

eq144.gif(2.88)

[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 le i le N, where N is fairly large, say,

eq145.gif(2.89)

are nearly the same;

[150]  (b) all components of dfi, 1 le i le N have vanishing mathematical expectations and the same variance; and

[151]  (c) all components dfi, 1 le i le N are uncorrelated with each other.

[152]  Fourth type of the a priori information is the knowledge of the functional

eq146.gif(2.90)

establishing the preference relation on the set of approximate solutions of SLAE (2.82); here, G is an Ntimes m matrix with m ll 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

eq147.gif(2.91)

the inequality

eq148.gif(2.92)

implies that the approximate solution x(2) is preferable.

[154]  In the first situation ( 1 < g2 le 2 ), the vector x must evidently satisfy the conditions

eq149.gif(2.93)

[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:

eq150.gif(2.94)

[157]  In the second approach, the vector

eq151.gif(2.95)

is first determined, and the vector x is then found as a sufficiently accurate approximate solution of the SLAE

eq152.gif(2.96)

As shown below, the second approach is preferable in the case of SLAE (2.82) of a large dimension (tentatively, min(N, M) ge5000 ).

[158]  The point is that the vector u can (and must!) be found from the solution of the constrained extremal problem

eq153.gif(2.97)

Furthermore, it is evident (and basically important!) that the conditions

eq154.gif(2.98)

can be replaced by the conditions

eq155.gif(2.99)

Therefore, constrained extremal problem (2.97) is replaced by the following one:

eq156.gif(2.100)

Using the method of Lagrangian multipliers, constrained extremal problem (2.97) corresponds to the family of unconstrained extremal problems

eq157.gif(2.101)

where

eq158.gif(2.102)

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

eq159.gif(2.103)

[160]  Since GG T is an Ntimes N matrix, it is appropriate to pass from SLAE (2.103) to the extended SLAE of second kind

eq160.gif(2.104)

where

eq161.gif(2.105)

eq162.gif(2.106)

and

eq163.gif(2.107)

Here EN is the unit Ntimes N -matrix.

[161]  The value l is always sufficiently large,

eq164.gif(2.108)

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

eq165.gif(2.109)

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:

eq166.gif(2.110)

[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

eq167.gif(2.111)

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

eq168.gif(2.112)

where

eq169.gif(2.113)

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

eq170.gif(2.114)

[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,

eq171.gif(2.115)

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

eq172.gif(2.116)

[171]  It is clear that the most rational approach is to use new iterative methods. The criterion of terminating the process of successive iterations

eq173.gif(2.117)

should be taken in the form

eq174.gif(2.118)

where Crit is an a priori given number related to the value

eq175.gif(2.119)

Namely, it is appropriate to take

eq176.gif(2.120)

[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:

eq177.gif(2.121)

where the quantity p2 is given and

eq178.gif(2.122)

[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

eq179.gif(2.123)

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

eq180.gif(2.124)

where

eq181.gif(2.125)

are the Lagrangian multipliers accounting for the equality conditions in problem (2.121).

[179]  The solution of problem (2.124) satisfies the SLAE

eq182.gif(2.126)

where

eq183.gif(2.127)

and, because G is an Ntimes m matrix, B is an Mtimes M matrix. Evidently, SLAE (2.126) should be transformed into the extended SLAE of second kind

eq184.gif(2.128)

where

eq185.gif(2.129)

eq186.gif(2.130)

and

eq187.gif(2.131)

Here, EN+m and EM are unit matrices of the sizes (N+m)times(N+m) and Mtimes M, respectively.

[180]  As in the case of SLAE (2.128), the vector xl, m can be represented in the form

eq188.gif(2.132)

where xl, 1 and xl, 0 are determined through the respective SLAE solutions

eq189.gif(2.133)

eq190.gif(2.134)

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

eq191.gif(2.135)

where Ll is the lower triangular matrix of the size (N + m + M)times(N + m + M) with unity diagonal elements and Rl is the upper triangular matrix of the same size. The determination of l and m values ensuring the validity of the equality conditions in (2.121) is quite similar to the case described in the preceding paragraph (in the case of the filtering problem of finding the vector uapprox f ).

[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).


RJES

Citation: Strakhov, V. N. (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 Sciences

Powered by TeXWeb (Win32, v.2.0).