CN113238284A - Gravity and magnetic fast forward modeling method - Google Patents

Abstract

The invention discloses a gravity and magnetic fast forward modeling method. First, an observation height and a geophysical model are determined, the geophysical model is discrete, a virtual line observation grid/virtual surface observation grid is proposed, and a core matrix containing multiple BTTB matrices is obtained, A series of kernel matrices corresponding to different layers are constructed by using multi-scale Haar wavelets, so as to realize the fast forward method of gravity and magnetism based on the multi-level method. Compared with the traditional multi-scale Haar wavelet to construct the kernel matrix of different layers, the original calculation and storage of the kernel matrix is required to cause the problem of excessive memory consumption. The memory consumption is reduced, and the calculation efficiency of gravity and magnetic forward calculation is improved.

Description

A Gravity-Magnetic Fast Forward Method

technical field

The invention belongs to the technical field of geophysics and exploration, in particular to the technical field of gravity and magnetic exploration.

Background technique

Fast and efficient forward calculation is the basis of large-scale potential field inversion. At present, the validity and timeliness of gravity and its gradient forward method are the important reasons restricting the inversion and interpretation of large-scale gravity data.

In order to solve the geophysical multi-scale forward modeling problem, the multi-grid method is generally introduced, which first defines the fine grid, and then defines the next-level coarse grid according to the fine grid. and high-frequency components) are assigned to grids of different thicknesses for solving. According to the different ways of defining transfer operators between grids, multigrid methods can be divided into geometric multigrid methods and algebraic multigrid methods. The geometric multigrid method does not require any preprocessing, and can accurately and efficiently complete the calculation for boundary value problems with uniform grids and uniform physical properties, but it is difficult to apply to unbounded open-domain problems, and its too sparse and coarse grids are difficult to truly reflect The distribution of physical properties, and the transfer operator between grid layers cannot accurately transfer information between different grid layers, so the multigrid algorithm usually loses its inherent efficiency; algebraic multigrid abandons the concept of geometric grids , a virtual grid defined entirely based on algebraic methods is introduced, which does not need to clarify the geometric and physical meanings of each layer of grids, and has the advantages of small storage, high convergence accuracy and low calculation time. The serial nature of the strategy" hinders its application to large-scale parallel computing.

Based on the above description, there is an urgent need for a Gravity-Magnetic fast forward method to solve the problem of excessive memory consumption caused by the traditional calculation and storage of kernel matrices for constructing different layers of kernel matrices.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a method for fast forward calculation of gravity and magnetism, by which it does not need to consume excessive memory when constructing different layers of kernel matrices, and realizes fast and efficient calculation of gravity and magnetic fast forward calculation.

In order to solve the above-mentioned technical problem, the concrete technical scheme of a kind of gravity and magnetic fast forward method of the present invention is as follows:

A Gravity-Magnetic Fast Forward Modeling Method Based on Multi-level Method, comprising the following steps:

S1, determine the multi-scale grid layer n, discrete geophysical model;

S2. Determine the calculation formula of gravity and magnetic potential field;

S3, set a virtual line observation grid to obtain a kernel coefficient matrix G _q,m ;

S4. Set the virtual surface observation grid to obtain the kernel coefficient matrix

S5. Use multi-scale Haar wavelets to construct kernel matrices corresponding to different layers, and realize fast forward modeling of gravity and magnetism based on multi-level method

Preferably, step S1 adopts the point element method to divide the geophysical model into a plurality of upright hexahedrons using parallel axial section planes.

Preferably, the calculation formula of the gravity and magnetic potential field in step S2 is in the form of a matrix:

d=Gm

In the formula, d is the observation data, all of which are vectors of length N _d ; m is the physical parameter vector of length N _m , that is, the physical parameters of the physical grid, such as density value, magnetization or magnetic susceptibility value, etc.; G is the corresponding Forward calculus operator, a matrix of size N _d ×N _m ; N _m and N _d are the number of grids and the number of observation data points respectively; N _m =n _x × _ny ×n _z , N _d =n _x × _ny , where n _x , _ny and _nz are the section fractions of the model along the three axes, respectively.

Preferably, the observation data d is g _z , g _xx , g _xy , g _xz , g _yy , g _yz or g _zz .

Preferably, in step S3, G _q,m is the kernel coefficient matrix of the q-th survey line and the corresponding near-surface m-th column grid, and the product of G _q,m and the corresponding physical property unit m _q,m G _{q, m} *m _q,m is:

和F为傅里叶变换对；作为替代核系数矩阵而需要存储的部分

为：In the formula,

and F are the Fourier transform pair; the part that needs to be stored as a substitute for the kernel coefficient matrix

for:

为第t层观测点及所对应的第n层物性网格观测方案的核系数矩阵，

与对应的物性单元

的乘积

为：Preferably, in step S4

is the kernel coefficient matrix of the t-th layer observation point and the corresponding n-th layer physical property grid observation scheme,

with the corresponding physical unit

product of

for:

The part that needs to be stored as a substitute for the kernel coefficient matrix

Preferably, the conversion of the adjacent two-layer kernel coefficient matrices is:

G ⁱ⁺¹ =R ⁱ G ⁱ P ⁱ (14)

In the formula, R is the restriction operator, and P is the transfer operator or the interpolation operator.

Preferably, in step S5, the multi-scale Haar wavelet is used to construct the kernel matrix corresponding to different layers as follows:

Choose Haar wavelet as R ^T and P, then:

In the formula, j, k=1, 2, W and W ^T are the Haar wavelet matrix forward transformation and inverse transformation, respectively, the superscript i represents the layer number of preprocessing, i=0 corresponds to the finest grid, and i=n is Thickest grid.

为

的Toeplitz向量为

Preferably, in step S5

for

The Toeplitz vector is

为W_j小波基。In the formula,

is W _j wavelet basis.

Beneficial effects:

The invention makes full use of the characteristics that the gravity and magnetic potential field forward kernel matrix is a special structure matrix and a multi-scale Haar wavelet, converts the large-scale gravity and magnetic forward kernel matrix to the wavelet domain, and realizes the Toeplitz-based wavelet multi-scale domain fast normalization. play. In addition, it has the following advantages:

1. A virtual line-line observation grid/virtual surface-surface observation grid is proposed, and a kernel matrix containing multiple BTTB matrices is obtained.

2. Propose a kernel matrix method that can construct different layers only by using Toeplitz vectors, to avoid the problem of excessive memory consumption caused by the traditional calculation and storage of kernel matrices required to construct different layers of kernel matrices.

3. Use Haar multi-scale wavelet to ensure that the kernel matrix of different layers still has a special structure, so as to obtain efficient computing performance, and then obtain higher-precision forward calculation results of gravity and magnetism, gravity and magnetism vector and gravity and magnetic gradient tensor.

Description of drawings

Fig. 1 is a flow chart of a method for optimizing gravity and magnetic fast forward modeling provided by the present embodiment;

2 is a schematic diagram of a virtual line observation scheme provided by the present embodiment;

FIG. 3 is a schematic diagram of a virtual surface observation scheme provided by the present embodiment;

FIG. 4 is a schematic diagram of a forward modeling model provided in this embodiment.

Detailed ways

In order to better understand the method and steps of the present invention, a method for fast forward inversion of gravity and magnetism of the present invention will be described in further detail below with reference to the accompanying drawings.

As shown in Figure 1 to Figure 4, the gravity and magnetic fast forward method provided by the present invention is based on a multi-level method, and its theory and derivation are as follows:

In the Cartesian coordinate system, there exists a three-dimensional body with residual mass m, volume V, and residual density ρ. The point element method is used to divide the geophysical model into a large number of upright hexahedrons using a series of parallel axonometric planes. There is no analytical singularity forward calculation formula for the total gravity gradient tensor of any observation point P to any cuboid [ξ ₁ →ξ ₂ , η ₁ →η ₂ , ζ ₁ →ζ ₂ ] in the underground space, taking the gravity g _z as an example :

In the formula,

μ _ijk =(-1) ^j+j+k ,

x _i =x-ζ _i , y _j =y-η _j , Z _k =z-ζ _k .

For the abnormal response of all observation points, the abnormal response of each upright hexahedron in the underground space to the corresponding observation point can be calculated one by one according to the superposition principle. Therefore, the forward calculation of gravity and gravity tensor can be written in matrix form:

d=Gm (1)

In the formula, d is the observation data, which can be specifically g _z , g _xx , g _xy , g _xz , g _yy , g _yz and g _zz , etc., all of which are vectors of length N _d ; m is a vector of physical property parameters of length N _m , that is, the physical parameters of the physical grid, such as density value, magnetization or magnetic susceptibility value, etc.; G is the corresponding forward operator, which is a matrix of size N _d × N _m ; N _m and N _d are the meshes, respectively and observed data points. N _m =n _x × _ny ×n _z , N _d =n _x × _ny , where n _x , _ny and _nz are the section numbers of the model along the three axes, respectively.

For traditional gravity and magnetic exploration, in the Cartesian coordinate system, taking the gravity field forward calculation as an example, when the calculation formula of the forward calculation kernel function is determined, the kernel function is only related to the spatial relative position of the observation point and the corner point of the grid cell. related. Here, the triaxial meshing fractional notation <*,*,*> is introduced, which is a function of n _x , n _y . The j-th physical property grid Q and the i-th observation point P are described by the three-axis division points along the Cartesian coordinate system as <l,m,n> and <p,q,t> respectively, then:

In the formula, <1,1,1> is the calculation origin.

It can be seen from equation (2) that, by using the symmetry of the kernel function memory, the calculation of the kernel function of any observation point in any grid can be calculated by converting to the calculation origin, which greatly simplifies the calculation of certain symmetric characteristics. components, but not all components have this property.

As shown in Figure 2, a virtual line observation scheme is set. The so-called line observation scheme here is the observation scheme of a single column of physical property grids in the geophysical model corresponding to a single survey line. Therefore, the equivalent geometric frame calculation formula based on translation equivalence is:

In the formula, 1≤l≤n _x , 1≤m≤n _x , _1≤p≤ny , 1≤q≤n _x , 1≤n≤n _z , t≥1.

For conventional exploration methods, observation data d is two-dimensional, that is, t=1; for isodimensional inversion, observation data d is three-dimensional, and t≥1. Using formula (3) will greatly simplify the forward calculation of the different components.

In the formula, i'=(q-1) _nx +(t-1)nxny, and j' ₌ (m-1) _ny +( _n - ₁ ) _nxny .

In Equation (3), p and l move from 1 to n _x along the x-axis with Δp and Δl (and Δp=Δl) as steps, respectively. From Equation (4), it can be known that:

According to formula (5), for any q and m, the diagonal line with equal elements can be found in the matrix on the left side of formula (4). For example, when q=1 and m=1, Δp goes from 1 to n _x along the line The x-axis starts to move, that is:

In a similar way, other diagonals can be constructed, that is, the kernel coefficient amplitudes of any diagonal are the same. Therefore, the matrix on the left side of equation (4) can be expressed by 2n _x -1 elements, which can be rewritten as:

In the formula, G _q,m corresponds to the qth survey line and the corresponding kernel coefficient matrix of the mth column of the grid near the surface, and a _-p and a _l are the diagonal lines in the lower/upper triangular matrix of formula (4). the first element of .

According to advanced mathematical knowledge, the product G _q,m *m _q,m of G _q,m and the corresponding physical property unit m _q ,m can be written as:

和F为傅里叶变换对；作为替代核系数矩阵而需要存储的部分

为：In the formula,

and F are the Fourier transform pair; the part that needs to be stored as a substitute for the kernel coefficient matrix

for:

Here, formula (8) is defined as a function T(G _{q,m ,} m _{q,m )} . As far as the sum of the products of a single vector v is concerned, by deriving similar formulas (4) to (8), it can be obtained:

∑ _n=1 G _q,m,n v=(∑ _n=1 G _q,m,n )v (9)

The above formula is mainly used in the calculation of G ^T d in the potential field inversion, which can further greatly speed up the implementation of the inversion.

As shown in Figure 3, a virtual surface observation scheme is set. The so-called surface observation plan here is the observation plan composed of the observation surface composed of multiple survey lines and the corresponding single-layer physical property grid in the geophysical model. This plan does not exist in real geophysical exploration and is only for derivation. based on the fast forward calculation method.

The relationship between other survey lines and each column of physical property grids is constructed by formula (7), and it is extended to the t-th layer observation point and the corresponding kernel coefficient matrix of the n-th layer physical property grid observation scheme

的计算公式：Similarly, the derivation process similar to formula (8) can be used to derive

Calculation formula:

The part that needs to be stored as a substitute for the kernel coefficient matrix

和

因G_q,m为Toeplitz矩阵，故而G_q,m和

也为Toeplitz矩阵。For traditional potential field exploration, there are three types of observation schemes: line observation scheme, surface observation scheme and iso-dimensional inversion scheme, whose kernel coefficient matrices are G _q,m ,

and

Since G _q,m is a Toeplitz matrix, G _q,m and

Also a Toeplitz matrix.

Aiming at the bottleneck problem that the point-element method is difficult to calculate or store the kernel matrix when dealing with large-scale, multi-area, and large-volume gravitational and magnetic gradients and their tensor data, many methods based on multi-scale pyramids and algebraic multi-grid are used for reference. The idea of scale layering is to process large-scale, multi-area and large-volume data step by step, so that traditional computers can also process large-scale, multi-area and large data volume data. The formula (1) can be rewritten as:

G ⁱ m ⁱ =d ⁱ , 0≤i≤n (13)

In the formula, the superscript i represents the layer number of the preprocessing, i=0 corresponds to the finest grid, that is, G ⁰ =G, and i=n is the coarsest grid.

The transformation of adjacent two-layer kernel coefficient matrices can be written as:

G ⁱ⁺¹ =R ⁱ G ⁱ P ⁱ (14)

In the formula, R is the restriction operator, P is the transfer operator, or the interpolation operator, that is, the transfer of the kernel function from the fine grid to the coarse grid, and vice versa.

When choosing Haar wavelet as R ^T and P:

Rewrite equation (13) to the frequency domain, then:

和

where W and W ^T are the Haar wavelet matrix forward transformation and inverse transformation, respectively,

and

Here, equation (15) is written in block matrix multiplication form using W ₁ and W ₂ :

和

in,

and

可以被证明为典型的Toeplitz矩阵和BTTB矩阵的组合,按虚拟线线观测方案将Gⁱ展开为矩阵形式：For j,k=1,2,

It can be shown to be a combination of a typical Toeplitz matrix and a BTTB matrix, expanding G ⁱ into matrix form according to the virtual line observation scheme:

的Toeplitz向量为：here the matrix

The Toeplitz vector is:

Furthermore, there are:

表示矩阵E和F的张量积，vec(E)为对矩阵E的向量化运算。in,

Represents the tensor product of matrices E and F, and vec(E) is a vectorized operation on matrix E.

因而在内存消耗上并没有改进。为此，将式(18)进行改写获得

的Toeplitz向量

Since the calculation formula (18) needs to use the original kernel matrix

So there is no improvement in memory consumption. To this end, formula (18) can be rewritten to obtain

Toeplitz vector

为W_j小波基。In the formula,

is W _j wavelet basis.

In order to illustrate the present invention in conjunction with the embodiment, a theoretical model is established, as shown in FIG. 1 , the specific implementation process of the present invention based on the multi-level method gravity and magnetic fast forward method is as follows:

1) Determine the number of layers of the multi-level method, such as n=3, divide the geophysical model according to the fact that _nx = _ny is divisible by ²³ ;

2) Determine the calculation formula of the gravity and magnetic potential field, where the gravity g _z is taken as an example;

3) with the virtual line observation scheme, calculate the Toeplitz vector t ⁰ that needs to be used to generate the kernel matrix;

4) Construct Toeplitz vectors ^ti at different layers;

5) Constructing multi-scale wavelet operators in different layers

6) In the wavelet domain, realize fast forward modeling of different block matrices

It can be understood that those skilled in the art know that various changes or equivalent substitutions can be made in this method without departing from the spirit and scope of the present invention. In addition, in accordance with the teachings of the present invention, the method may be modified to adapt a particular situation without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited by the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of the present application fall within the protection scope of the present invention.

Claims (9)

1. a fast forward method for gravity and magnetism, is characterized in that, comprises the following steps: S1, determine the multi-scale grid layer n, discrete geophysical model; S2. Determine the calculation formula of gravity and magnetic potential field; S3, set a virtual line observation grid to obtain a kernel coefficient matrix G _{q, m} ; S4. Set the virtual surface observation grid to obtain the kernel coefficient matrix

S5. Use multi-scale Haar wavelets to construct kernel matrices corresponding to different layers, and realize fast forward modeling of gravity and magnetism based on multi-level method

2 . The method according to claim 1 , wherein in step S1 , the point element method is used to divide the geophysical model into a plurality of upright hexahedrons using parallel axis-section planes. 3 . 3. a kind of gravity and magnetic fast forward method according to claim 1, is characterized in that, in step S2, the calculation formula of gravity and magnetic potential field is matrix form: d=Gm In the formula, d is the observation data, all of which are vectors of length N _d ; m is the physical parameter vector of length N _m , that is, the physical parameters of the physical grid, such as density value, magnetization or magnetic susceptibility value, etc.; G is the corresponding Forward calculus operator, a matrix of size N _d ×N _m ; N _m and N _d are the number of grids and the number of observation data points respectively; N _m =n _x × _ny ×n _z , N _d =n _x × _ny , where n _x , _ny and _nz are the section fractions of the model along the three axes, respectively. 4 . The method according to claim 3 , wherein the observation data d is g _z , g _xx , g _xy , g _xz , g _yy , g _yz or g _zz . 5 . 5. a kind of gravity-magnetic fast forward method according to claim 3, is characterized in that, in step S3, G _{q, m} is the kernel coefficient matrix of the qth survey line and the corresponding near-surface mth column grid , G _{q, m} and the corresponding physical property unit m _{q, m} product G _{q, m} *m _{q, m} is:

In the formula,

and F are the Fourier transform pair; the part that needs to be stored as a substitute for the kernel coefficient matrix

for:

6. a kind of gravity-magnetic fast forward method according to claim 5, is characterized in that, in step S4

is the kernel coefficient matrix of the t-th layer observation point and the corresponding n-th layer physical property grid observation scheme,

with the corresponding physical unit

product of

for:

The part that needs to be stored as a substitute for the kernel coefficient matrix

7. a kind of gravity-magnetic fast forward method according to claim 6, is characterized in that, the conversion of adjacent two layers of kernel coefficient matrix is: G ⁱ⁺¹ =R ⁱ G ⁱ P ⁱ (14) In the formula, R is the restriction operator, and P is the transfer operator or the interpolation operator. 8. a kind of gravity-magnetic fast forward method according to claim 7 is characterized in that, in step S5, using multi-scale Haar wavelet to construct the kernel matrix corresponding to different layers is: Choose Haar wavelet as R ^T and P, then:

In the formula, j, k=1, 2, W and W ^T are the Haar wavelet matrix forward transformation and inverse transformation, respectively, the superscript i represents the layer number of preprocessing, i=0 corresponds to the finest grid, and i=n is Thickest grid. 9. a kind of gravity-magnetic fast forward method according to claim 8, is characterized in that, in step S5

for

The Toeplitz vector is

In the formula,

is W _j wavelet basis.

Images