Device For UnderGround Exploration

 

Due to its low cost and convenient construction method, gravity prospecting is favored

Widely used in large-scale geological anomaly exploration and large-scale ore prospecting

Census, and small-scale density 3D geological modeling work.

There are two commonly used inversion methods at present, 2.5D combined with 3D interface inversion.

[1-2] and 3D physical property inversion. By performing three-dimensional

Grid division, inversion of the physical properties of each grid unit, so as to obtain geological

Physical property distribution and morphological characteristics of abnormal body[3]

.

For the 3D inversion of targets with complex structures, due to the

According to the limitations of calculation and storage capacity, three-dimensional physical property reflection is mainly used.

play. Even so, 3D forward simulation and inversion of gravity field calculations

still faces a large amount of data processing, using the traditional linear

Inversion, often the calculation of millions of grid points will exceed the personal computer

The computer's memory processing range. Due to computer memory and CPU

The bottleneck of computing efficiency greatly limits the data processing and interpretation work.

quality of work. In order to solve this problem, in recent ten years, parallel computing

It has become a mainstream research direction. Chen et al. [4-5] adopted a method based on

The CUDA code of C language realizes the three-dimensional forward imaging of gravity;

Lin Wei et al[6] used GPU genetic algorithm to realize the three-dimensional ocean gravity

Inversion; Cuma et al. [7-8] realized the CPU-based OpenMP algorithm

Acceleration with GPU-based OpenACC algorithm acceleration, providing in-

The realization of gravity inversion calculation on high-performance computer clusters such as servers

method. From the hardware structure point of view, compared with only a few, several

CPUs with ten cores, GPUs with hundreds or even thousands of cores

Simple operations on large amounts of data, such as gravity inversion large matrix operations,

Has a natural advantage.

The author uses the traditional Fortran language on the PGI Fortran platform

language and CUDA language mixed, the reweighted regularization conjugate gradient

(Re-weighting Regularized Conjugate Gradient,

RRCG) 3D gravity field inversion algorithm is reorganized so that it can be used in

GPUs on personal desktop and notebook computers

Parallel computing, the use of regularization makes the inversion results stable and reliable, with

physical meaning; abandoning the storage sensitivity matrix A, thus making the computational efficiency

The efficiency and data processing volume are improved at the same time; multiple thread calls are introduced

By searching with the sensitivity range (footprint), the storage expansion and

Inversion accelerated. In order to verify the correctness of the program, we use the same

The MATLAB traditional code of the algorithm was compared and tested. the party

The method can effectively solve the reflection of hundreds of millions of grid points on a personal computer

to solve the problem of 3D gravity inversion with a large amount of data

A new way of thinking.

1 Forward modeling and inversion

1.1 Forward Calculation Algorithm

Although the analytical solution for the volume integral of the gravitational component.Discretized single-point integration of : Divide the subsurface medium into uniform cubic

volume grid, and assuming that the internal density of the grid is uniform and uniform, after grid separation After bulking, the commonly used vertical component calculation formula in gravity measurement discrete form. Among them: γ represents the gravitational constant; Nm represents the constant total number of square grids; k represents the number of grid points; vector r′ = ( x′, y′, z′), r = (x, y, z) represent the observation point and the source grid center respectively

The position of the point; z and z′ represent the vertical distance between the observation point and the source center point respectively. Δx, Δy, Δz respectively denote the partitioned cube network

The length of the lattice along the three directions of the relative coordinate axis; ρ represents the cube network grid density.

Under this condition, the process of calculating the gravity anomaly value of the observation point can be it can be regarded as the relationship between a linear operator A and the density distribution vector operation

d = Am, (2)

Among them: d is expressed as a spatial distribution of observation points with a length of Ndata vector, m is expressed as a space density distribution vector with length Mcell Quantity, A represents the linear operator of the vertical component of gravity, in such a linear Under the problem, A is also the sensitivity matrix in the inversion problem.

1.2 Inversion algorithm

Considering the speed of parallel computing, the author uses the

Reweighted regularized conjugate gradient (RRCG) algorithm, Tikhonov

The introduction of regularization can suppress the problems caused by data noise and multi-solution

The resulting inversion results are unstable and the physical meaning is not clear[9-10]

instead of

Different weighting factors can be between smooth solutions and solutions with sharp boundaries

Balance, reweighting (Re-weighting) can be in the iterative process

Dynamically adjust the range of the regularization factor α, which can guarantee the quality of the data

Accelerate the convergence speed under the premise of quality. The author only introduces and parallel computing

Related algorithm flow, for the selection of weighting factor and regularization factor

The selection and update methods will not be introduced anymore, please refer to the literature [9] for details.

.

(7)

Among them: rn represents the residual error vector; n represents the number of iterations; l

n table

is the iterative direction vector of the gradient method, and ln represents the iterative direction vector of the conjugate gradient method

vector; F represents the Frechet derivative matrix; in the linear inversion problem

where F=A; kn represents the iteration step size of each iteration; αn is the

Regularization parameters for step iterations.

In the traditional serial RRCG inversion, calculation and storage sensitive

The degree matrix A is read repeatedly, and each iterative step is calculated sequentially.

Stop until the iteration meets the number condition or error condition. A moment

The matrix is a large matrix related to the multiplication of data and model scales.

Time-consuming and difficult to store and read in large-scale 3D inversions, becoming a limitation

important factor in conventional RRCG linear inversion.

2 GPU Algorithm Design

2.1 GPU and CUDA

GPU (Graphics Processing Unit), that is, graphics processing unit

Yuan, the computing method is highly parallel, and the computing power is powerful, which was adopted by the early science

It is used by researchers to improve the calculation speed of general computing. And nearly ten years

Later, due to the emergence of CUDA, the GPU has become an unnecessary

A general-purpose parallel computing tool that can be used to master the basics of graphics computing

with [12]

. CUDA (Compute Unified Device Architecture)

It is a general-purpose GPU programming language framework provided by NVIDIA, currently

It can support the use of multiple language platforms such as C, Fortran, and Matlab.

Use the CUDA statement to call the GPU computing unit for parallel computing.

CUDA is the most suitable for the linear 3D inversion of the potential field.

One of the means of parallel computing. First, using the RRCG inversion method,

Its process is an iterative algorithm, and each calculation step can be simplified as

Perform simple operations on independent matrices and vectors for parallel computing

provided the conditions; secondly, compared with the traditional CPU serial or parallel,

The number of computing cores of the GPU has an absolute leading advantage, compared with the target

The 16 or 32 computing cores of former high-end personal computers,

Personal graphics cards can have 400-600, or even more than 2000

GPU computing cores; finally, compared to another recently released

GPU computing standard OpenACC, although CUDA requires more