Method for automatically constructing meshes for model inversion using finite element methods for magnetotelluric data

Abstract

FIELD: geophysics.SUBSTANCE: methods of processing and interpreting magnetotelluric sounding data. The method comprises the following stages: using measured magnetotelluric data obtained at a plurality of observation pickets, the parameters of electromagnetic field penetration into the geoelectric medium are determined for different probing frequencies, on the basis of which the minimum thickness of the near-surface layer and the maximum depth of the computational domain are automatically set; for each observation picket, an individual vertical geometry of the layers of the computational domain is formed, based on their totality, the distribution of the density of the partitioning of the computational domain by depth is determined and on this basis a single vertical adaptive mesh of the model partitioning is formed; the spatial position of the upper boundaries of the finite elements is specified taking into account the relief and their smooth transition to the horizontal position at a given depth is ensured; the horizontal partitioning of the computational domain is formed adaptively with a local change in the sizes of the finite elements depending on the differences in the vertical geometry of the layers, strongly vertically elongated cells of the central part of the model are combined in pairs to increase the degree of isotropy of the cells at depth, the formed adaptive structured mesh is used for the numerical solution of the inverse problem of magnetotelluric sounding by the finite element method.EFFECT: increase in the accuracy of approximation of the electromagnetic field in a geoelectric medium and an increase in the stability of the solution of the inverse problem of magnetotelluric sounding.8 cl

Description

[0001] FIELD OF TECHNOLOGY

[0002] The claimed technical solution generally relates to the field of geophysics, namely to methods of processing and interpreting magnetotelluric sounding data.

[0003] In particular, the claimed solution relates to a method for constructing computational grids used in the numerical solution of inverse problems of magnetotelluric sounding using the finite element method.

[0004] STATE OF THE ART

[0005] When solving inverse problems of magnetotelluric sounding, it is necessary to construct a computational grid that ensures a correct approximation of the electromagnetic field in the geoelectric medium.

[0006] Various methods for constructing computational grids are known from the prior art.

[0007] There are known methods for constructing computational grids in which the computational domain is divided into a grid with a constant horizontal step and a logarithmically increasing vertical step.

[0008] This approach is described, for example, in the source:

[0009] Constable SC, Parker RL, Constable CG Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 1987, Vol. 52, No. 3, pp. 289-300.

[0010] The disadvantages of this approach are:

[0011] lack of grid adaptation to the structure of measured data;

[0012] insufficient mesh detail near observation points;

[0013] Sharp increase in the number of elements when using a fine uniform grid.

[0014] A known method of grid construction involves determining the layer depths based on the penetration depth of the electromagnetic field. The penetration depth is determined through the skin layer thickness, which depends on the probing frequency and the electrical resistance of the medium. This approach is described, for example, in the following sources (Berdichevsky M.N., Dmitriev V.I. Magnetotelluric Sounding of Horizontally Homogeneous Media. Moscow: Nedra, 1991) and Ulugergerli E.U., Candansayar M.E. Automated mesh design for 2-D magnetotelluric inversion. Geophysical Journal International, 2002, Vol. 148, pp. 403-414).

[0015] The disadvantages of this method are:

[0016] grid adaptation is performed only vertically;

[0017] The horizontal variability of the geoelectric structure is not sufficiently taken into account;

[0018] When there are 2D or 3D inhomogeneities, the simulation accuracy decreases.

[0019] In addition, methods for constructing grids are known in which the computational grid is condensed in the areas where magnetotelluric observation points are located.

[0020] A similar approach is described, for example, in the work: Tietze K., Ritter O., Egbert G.D. 3-D joint inversion of the magnetotelluric phase tensor and vertical magnetic transfer functions. Geophysical Journal International, 2015, Vol. 203, pp. 1128-1148.

[0021] The disadvantages of this method are:

[0022] The cell size is determined only by the configuration of the observation network;

[0023] The grid may be over-detailed in areas with simple geoelectric structure;

[0024] Inside conductive objects, the detail may be insufficient.

[0025] The technical solution proposed for consideration is aimed at eliminating the shortcomings of the current level of technology.

[0026] SUMMARY OF THE INVENTION

[0027] The technical problem, which the declared technical solution is aimed at solving, is the creation of a new method for automatically constructing meshes for dividing a model for inversion of magnetotelluric data using finite element methods.

[0028] The technical problem that the claimed solution is aimed at solving is to increase the accuracy and stability of the solution of the inverse problem of magnetotelluric sounding while simultaneously reducing the computational costs when constructing the computational grid.

[0029] The technical result consists in increasing the accuracy of approximation of the electromagnetic field in the geoelectric medium, increasing the stability of the solution of the inverse problem of magnetotelluric sounding.

[0030] An additional technical result is the automatic determination of the computational grid parameters.

[0031] The specified technical result is achieved through a method for automatically constructing meshes for dividing a model for inverting magnetotelluric data using finite element methods, which contains the following stages:

[0032] based on measured magnetotelluric data obtained at multiple observation points, physically determined parameters of the penetration of the electromagnetic field into the geoelectric environment are determined for different sounding frequencies;

[0033] Based on the specified parameters, the minimum thickness of the surface layer and the maximum depth of the calculation area are automatically set;

[0034] For each observation point, an individual vertical geometry of the layers of the calculation area is formed, corresponding to the effective penetration depths of the electromagnetic field;

[0035] based on the set of individual vertical geometries of the layers, the distribution of the density of the partitioning of the calculation area by depth is determined and, on this basis, a single vertical adaptive mesh of the model partitioning is formed;

[0036] If there is a relief on the earth's surface, the spatial position of the upper boundaries of the finite elements is specified taking into account the relief and their smooth transition to a horizontal position at a given depth is ensured;

[0037] The horizontal partitioning of the computational domain is also formed adaptively with a local change in the sizes of finite elements depending on the differences in the vertical geometry of the layers corresponding to the different components of the electromagnetic field,

[0038] The strongly vertically elongated cells of the central part of the model are combined in pairs to increase the degree of isotropy of the cells at depth,

[0039] The resulting adaptive structured grid is used to numerically solve the inverse problem of magnetotelluric sounding using the finite element method.

[0040] In a particular embodiment, the minimum thickness of the near-surface layer is set proportional to the minimum value of the skin layer thickness determined for all pickets and frequencies.

[0041] In another particular embodiment, the maximum modeling depth is set proportional to the maximum value of the skin layer thickness determined for all pickets and frequencies.

[0042] In another particular embodiment of the implementation, the powers of the layers of a single partition grid are set based on the minimum powers of the adjacent layers among all local grids at a given depth.

[0043] In another particular embodiment, the level of grid flattening is set at a depth determined based on the difference in elevation of the earth's surface relief.

[0044] In another particular embodiment, horizontal mesh refinement is performed only in areas of densification, within which the upper layers of individual vertical meshes obtained from different components of the impedance tensor differ significantly.

[0045] In another particular embodiment, the horizontal step of mesh detailing in the areas of thickening is determined on the basis of the minimum thickness of the upper layers of individual vertical meshes obtained from different components of the impedance tensor.

[0046] In another particular embodiment, cells characterized by a vertical to horizontal size ratio greater than 2 are subject to horizontal merging.

[0047] DETAILED DESCRIPTION OF THE INVENTION

[0048] The following detailed description of the invention includes numerous implementation details to provide a clear understanding of the present invention. However, one skilled in the art will readily understand how the present invention may be used with or without these implementation details. In other instances, well-known methods, procedures, and components have not been described in detail to avoid obscuring the features of the present invention.

[0049] Furthermore, it will be clear from the foregoing description that the invention is not limited to the embodiment described. Numerous possible modifications, changes, variations, and substitutions, while preserving the spirit and form of the present invention, will be apparent to those skilled in the art.

[0050] The claimed technical solution proposes a method for automatically constructing meshes for model partitioning for inversion of magnetotelluric data using finite element methods.

[0051] Below is the calculation of the skin layer thickness.

[0052] For all experimental apparent resistivity curves obtained at each observation site, the skin depth δ is calculated at each frequency f for the principal components of the impedance tensor (Zxy, Zyx). The skin depth characterizes the penetration depth of the electromagnetic field into the medium and determines the vertical resolution of the method at a given frequency. The calculation is performed using the standard formula (Electroelectric Prospecting: A Handbook of Electrical Prospecting Practice for Students of Geophysical Specialties. Volume 1 - Tver: Polipress, 2018):

[0053] δ=500√(ρ / f), where

[0054] δ - skin layer thickness (m),

[0055] ρ - apparent electrical resistance (Ohm⋅m).

[0056] The calculation is performed for the impedance tensor components:

[0057] Zxy

[0058] Zyx.

[0059] Defining the basic parameters of the grid

[0060] Based on the distribution of skin layer thicknesses, the following is determined:

[0061] minimum power of the first layer h1=min(δ) / 10, which ensures sufficient detail in the zone of maximum sensitivity of the MT field.

[0062] In this case, the coefficient 1 / 10 can take values ​​in the range: 1 / 5 - 1 / 20.

[0063] The maximum simulation depth is defined as Zmax=max(δ), which ensures correct signal attenuation at large depths.

[0064] Building individual vertical grids

[0065] For each observation point, a vertical grid is formed in which the layer boundaries correspond to the effective penetration depths of the field

[0066] z_eff=δ(f) / √2.

[0067] Vertical grids are built separately for components:

[0068] Zxy

[0069] Zyx.

[0070] Determination of partition density

[0071] For each vertical grid, a local partitioning density is calculated, which characterizes the number of layers per decade of depth.

[0072] The local density for the k-th boundary is determined by the formula n(z k )=21g(z k+1 / z k-1 ).

[0073] The obtained values ​​are interpolated to produce a smooth depth function n(z) at each station. The final density curve is then drawn using the maximum values ​​of all these curves. The height of each cell in the final model should not differ from the height of the cell above by less than 1.05 and more than 1.20 times.

[0074] If there is a relief on the earth's surface, the following operations are performed:

[0075] 1. The minimum and maximum elevation marks of the surface are determined;

[0076] 2. The difference in relief height Δz_rel is calculated;

[0077] 3. Sets the grid flattening level.

[0078] The depth of flattening is defined as

[0079] z_executed = z_min - Δz_rel.

[0080] Below this depth, the layer boundaries are assumed to be horizontal.

[0081] The horizontal mesh partition is determined based on the analysis of the differences between the vertical meshes constructed for the Zxy and Zyx components.

[0082] If the relative difference in the thickness of the corresponding layers exceeds a coefficient of 1.5, local condensation of the horizontal grid is performed.

[0083] The horizontal grid step in the thickening zone is defined as:

[0084] h=α⋅h_min,

[0085] where α takes values ​​in the range of 0.5-2,

[0086] h_min - minimum thickness of the upper layers of vertical meshes.

[0087] At the bottom of the model, the most vertically elongated cells with an anisotropy coefficient (i.e., the ratio of the vertical to horizontal cell size) greater than 2 are combined into pairs. This procedure can be repeated several times, provided it does not result in an excessively abrupt (more than threefold) change in the horizontal dimensions of adjacent cells.

[0088] The claimed solution is implemented using a computing system that processes magnetotellurgical measurements and automatically constructs a computational grid used in numerical modeling of the electromagnetic field when solving the inverse problem of magnetotelluric sounding using the finite element method.

[0089] These application materials present a preferred implementation of the claimed technical solution, which should not be used as limiting other particular embodiments of its implementation that do not go beyond the scope of the requested scope of legal protection and are obvious to specialists in the relevant field of technology.

Claims

1. A method for automatically constructing meshes for model partitioning for inversion of magnetotelluric data using finite element methods, containing the stages at which: Based on the measured magnetotelluric data obtained at multiple observation points, physically determined parameters of the penetration of the electromagnetic field into the geoelectric environment are determined for different probing frequencies; based on the specified parameters, the minimum thickness of the surface layer and the maximum depth of the calculation area are automatically set; for each observation point, an individual vertical geometry of the layers of the calculation area is formed, corresponding to the effective penetration depths of the electromagnetic field; based on the set of individual vertical geometries of the layers, the distribution of the density of the division of the calculation area by depth is determined and, on this basis, a single vertical adaptive grid of the model division is formed; in the presence of a relief of the earth's surface, the spatial position of the upper boundaries of the finite elements is specified taking into account the relief and their smooth transition to a horizontal position at a given depth is ensured; the horizontal partitioning of the computational domain is also formed adaptively with a local change in the sizes of the finite elements depending on the differences in the vertical geometry of the layers corresponding to the different components of the electromagnetic field, the strongly vertically elongated cells of the central part of the model are combined in pairs to increase the degree of isotropy of the cells at depth, The resulting adaptive structured grid is used to numerically solve the inverse problem of magnetotelluric sounding using the finite element method.

2. The method according to paragraph 1, characterized in that the minimum thickness of the near-surface layer is set proportional to the minimum value of the skin layer thickness, determined for all pickets and frequencies.

3. The method according to paragraph 1, characterized in that the maximum modeling depth is set proportional to the maximum value of the skin layer thickness, determined for all pickets and frequencies.

4. The method according to paragraph 1, characterized in that the powers of the layers of a single partition grid are set on the basis of the minimum powers of adjacent layers among all local grids at a given depth.

5. The method according to paragraph 1, characterized in that the level of grid flattening is set at a depth determined on the basis of the difference in elevation of the earth’s surface relief.

6. The method according to paragraph 1, characterized in that the horizontal detailing of the grid is carried out only in the zones of densification, within which the upper layers of individual vertical grids obtained from different components of the impedance tensor differ significantly.

7. The method according to paragraph 1, characterized in that the horizontal step of grid detailing in the areas of densification is determined on the basis of the minimum thickness of the upper layers of individual vertical grids obtained from different components of the impedance tensor.

8. The method according to paragraph 1, characterized in that the cells subject to horizontal merging are those characterized by a ratio of the vertical size to the horizontal size greater than 2.