Fast inversion method and system for semi-airborne transient electromagnetic based on multi-gpu parallelism

Abstract

The application provides a semi-airborne transient electromagnetic rapid inversion method and system based on multi-GPU parallelism, and relates to the technical field of semi-airborne transient electromagnetic inversion. Necessary initial information and an initial model are used to determine the calculation task of a Jacobian matrix; the calculation task of the Jacobian matrix is split into multiple subtasks on a CPU memory, the multiple subtasks are transmitted to corresponding GPUs respectively, and the calculation of the subtasks is performed in parallel on the GPUs; the results of all the subtasks are spliced into a complete Jacobian matrix on the CPU memory; an updated model is obtained based on the complete Jacobian matrix; forward calculation is performed on the updated model until the root mean square error meets a threshold value, and the final model obtained is a three-dimensional inversion result. The application splits the calculation task of the Jacobian matrix into multiple subtasks on the CPU memory, the multiple subtasks are distributed to multiple corresponding GPU cards, and the process of solving the inversion model is accelerated through a multi-GPU parallel calculation scheme.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of semi-airborne transient electromagnetic three-dimensional inversion in geophysical exploration, and particularly relates to a semi-airborne transient electromagnetic fast inversion method and system based on multi-GPU parallel. BACKGROUND

[0002] The semi-airborne transient electromagnetic method (SATEM) uses a ground wire source or a loop source as an excitation source to generate a stable electromagnetic field (also known as a primary field) in space during the supply of a stable current. After the current is quickly turned off, different media with different electrical properties in space will generate corresponding induced fields (also known as secondary fields). The secondary field signals carrying medium information can be obtained by receiving coils, and then the spatial distribution of the underground electrical structure can be obtained through interpretation means.

[0003] The structured grid-based finite difference time domain (FDTD) method has high calculation speed and has been widely used to solve SATEM problems. The implicit FDTD algorithm has unconditional stability, so the time step distance of this method can exceed the maximum time step required by the Courant condition. For a calculation model with small grid size, a larger time step can be used in each iteration to improve the calculation efficiency. The back-Euler direct-splitting FDTD (BEDS-FDTD) method establishes a full-discrete equation of the electromagnetic field containing a three-diagonal linear system by introducing an auxiliary field.

[0004] The three-dimensional inversion method is to interpret the SATEM observation data by optimization method. The Gauss-Newton (GN) method commonly used in optimization algorithms has a second-order convergence speed, so it is widely used in SATEM three-dimensional inversion calculation. The GN method requires explicit calculation of the Hessian matrix to solve the model increment, which requires explicit calculation of the Jacobian matrix, and the calculation speed of the Jacobian matrix directly affects the efficiency of three-dimensional inversion. In the process of calculating the Jacobian matrix using the adjoint method, it is also necessary to solve multiple sets of linear equations. The characteristic of this linear system is that it shares the same coefficient matrix, and the adjoint field values are solved through multiple right-hand side terms.

[0005] The traditional Jacobian matrix calculation method solves a linear system by calling a library function, and has a slow calculation speed. According to the characteristics of the BEDS-FDTD algorithm, the adjoint forward equation based on the forward theory is still a three-diagonal linear system, and therefore also has parallel solving characteristics. In the prior art, there is still a lack of a calculation method for quickly solving the Jacobian matrix developed for the BEDS-FDTD fast forward algorithm. SUMMARY

[0006] To overcome the deficiencies of the prior art, the present application provides a multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion method and system, which splits the Jacobian matrix calculation task into multiple sub-tasks on the CPU memory, assigns the multiple sub-tasks to multiple corresponding GPU graphics cards, and accelerates the process of solving the inversion model through a multi-GPU parallel computing scheme.

[0007] To achieve the above object, one or more embodiments of the present application provide the following technical solutions:

[0008] The present application provides a multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion method in a first aspect.

[0009] The multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion method comprises the following steps:

[0010] Obtaining necessary initial information for semi-airborne transient electromagnetic inversion;

[0011] Establishing an initial model and determining a regularization factor of the initial model;

[0012] Determining the calculation task of the Jacobian matrix based on the necessary initial information and the initial model;

[0013] Splitting the Jacobian matrix calculation task into multiple sub-tasks on the CPU memory, transmitting the multiple sub-tasks to corresponding GPUs respectively, and performing parallel computation of the sub-tasks on the GPUs;

[0014] After the GPUs perform the respective sub-tasks, transmitting the calculation results back to the CPU memory, and performing computation of other sub-tasks until the execution of all sub-tasks is completed;

[0015] Splicing the results of all sub-tasks into a complete Jacobian matrix on the CPU memory according to the splitting order;

[0016] Obtaining an updated model based on the complete Jacobian matrix;

[0017] Performing forward calculation on the updated model to obtain a root mean square error, and determining whether the root mean square error meets a threshold value, if not, modifying the regularization factor and repeating the calculation of the Jacobian matrix until the root mean square error meets the threshold value, and the final model obtained is the three-dimensional inversion result.

[0018] The second aspect of the present application provides a multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion system.

[0019] The multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion system comprises:

[0020] A data acquisition module configured to acquire necessary initial information for semi-airborne transient electromagnetic inversion;

[0021] A model establishment module configured to establish an initial model and determine a regularization factor of the initial model;

[0022] A Jacobian matrix determination module configured to determine a calculation task of a Jacobian matrix based on the necessary initial information and the initial model;

[0023] A task allocation module configured to split the calculation task of the Jacobian matrix into a plurality of subtasks on a CPU memory, transmit the plurality of subtasks to corresponding GPUs respectively, and perform parallel calculation of the subtasks on the GPUs;

[0024] A parallel execution module configured to, after the GPUs perform the respective subtasks in parallel, transmit the calculation results back to the CPU memory, and perform calculation of other subtasks until the execution of all subtasks is completed;

[0025] A return splicing module configured to splice the results of all subtasks into a complete Jacobian matrix in a splitting order on the CPU memory;

[0026] A model updating module configured to obtain an updated model based on the complete Jacobian matrix;

[0027] A judgment loop module configured to perform forward calculation on the updated model to obtain a root mean square error, judge whether the root mean square error meets a threshold value, modify the regularization factor and repeatedly calculate the Jacobian matrix if the root mean square error does not meet the threshold value, until the root mean square error meets the threshold value, and the final model obtained is the three-dimensional inversion result.

[0028] The above one or more technical solutions have the following beneficial effects:

[0029] The application provides a multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion method and system, aiming at the characteristics of GN optimization method, task scheduling and load balancing of multiple CPU cores are realized through OpenMP, efficient parallel calculation of the GPU end is realized by combining OpenACC technology, and a heterogeneous parallel computing architecture of CPU master and multiple GPU cooperation is constructed.

[0030] According to the characteristics of the adjoint forward equation based on the BEDS-FDTD forward algorithm, the coefficient matrix corresponding to multiple right-hand sides (RHS) has the same structural characteristics, and the traditional solving process is improved to a batch processing mode.

[0031] According to the characteristics of the adjoint forward equation based on the BEDS-FDTD forward algorithm, the coefficient matrix corresponding to multiple right-hand sides (RHS) has the same structural characteristics, and the traditional solving process is improved to a batch processing mode.

[0032] Advantages of the additional aspects of the application will be partially given in the following description, partially will become obvious from the following description, or will be known by the practice of the application. BRIEF DESCRIPTION OF DRAWINGS

[0033] The drawings constituting a part of the specification of the application are used to provide further understanding of the application, the illustrative embodiments of the application and the description thereof are used to explain the application, and do not constitute improper limitation of the application.

[0034] Figure 1 The method flowchart of example one.

[0035] Figure 2 The data structure diagram of the adjoint forward equation of multiple RHS of example one.

[0036] Figure 3 Flow chart of Thomas algorithm for solving three-diagonal matrix of example one.

[0037] Figure 4 Acceleration effect diagram of ThomasBatch algorithm of example one relative to cuThomasBatch algorithm. DETAILED DESCRIPTION

[0038] It should be noted that the following detailed description is exemplary in nature and is intended to provide further description of the application. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs.

[0039] It should be noted that the terms used herein are only for the purpose of describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application.

[0040] In the case of no conflict, the embodiments in the present application and the features in the embodiments can be combined with each other.

[0041] Example one

[0042] The present embodiment discloses a multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion method, which is suitable for data interpretation of semi-airborne transient electromagnetic (SATEM). The initial model, the updated model and the final model in three-dimensional inversion calculation respectively represent the grid division model containing the preset before calculation, the calculated electrical parameter distribution in the inversion process and the finally obtained electrical parameter distribution. These models are used as carriers to carry the electromagnetic field diffusion. By continuously updating the electrical parameters in the model, the electromagnetic response characteristics at the observation points are obtained. Therefore, these models need to consider the following parameters: size of the transmitting source, current in the transmitting source, waveform of the transmitting source, number and size of the grid, distribution of the observation points, etc.

[0043] As shown in Figure 1 The multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion method comprises the following steps:

[0044] Obtaining necessary initial information for semi-airborne transient electromagnetic inversion;

[0045] Establishing an initial model and determining a regularization factor of the initial model;

[0046] Based on the necessary initial information and the initial model, determining the calculation task of the Jacobian matrix;

[0047] Splitting the Jacobian matrix calculation task on the CPU memory into multiple sub-tasks, transmitting the multiple sub-tasks to corresponding GPUs respectively, and performing parallel computation of the sub-tasks on the GPUs;

[0048] After the GPUs perform the respective sub-tasks in parallel, the computation results are transmitted back to the CPU memory, and other sub-tasks are executed for computation until all the sub-tasks are executed;

[0049] After the GPUs perform the respective sub-tasks in parallel, the computation results are transmitted back to the CPU memory, and other sub-tasks are executed for computation until all the sub-tasks are executed;

[0050] Based on the complete Jacobian matrix, an updated model is obtained;

[0051] The updated model is subjected to forward calculation to obtain a root mean square error, and it is determined whether the root mean square error satisfies a threshold value; if not, a regularization factor is modified and the Jacobian matrix is repeatedly calculated until the root mean square error satisfies the threshold value, and the final model obtained is the three-dimensional inversion result.

[0052] For example, the initial model can be set as a half-space model with uniform electrical medium distribution.

[0053] Further, the necessary initial information for the semi-airborne transient electromagnetic inversion includes the conductivity distribution of the initial model, the position and size of the transmitting source, the position and time point of the observation data, the grid size, the boundary condition, and the grid division rule.

[0054] Further, the initial model represents a grid division model containing a pre-set electrical parameter distribution before calculation, the updated model represents a grid division model containing an electrical parameter distribution calculated in the inversion process, and the final model represents a grid division model containing a finally obtained electrical parameter distribution.

[0055] Further, based on the necessary initial information and the initial model, the calculation task of the Jacobian matrix is determined, specifically:

[0056] The Jacobian matrix is represented as:

[0057]

[0058] In the formula, d represents the observation data, N r and N m respectively represent the number of observation data and the number of inversion parameters; represents the i-row j-column Jacobian matrix; d i represents the i-th observation data; m jrepresents the conductivity of the i-th inversion grid.

[0059] Further, the Jacobian matrix calculation task on the CPU memory is split into multiple sub-tasks, and the multiple sub-tasks are transmitted to corresponding GPUs respectively, and the calculation of the sub-tasks is performed in parallel on the GPUs, specifically including:

[0060] The Jacobian matrix is expressed in the form of a block matrix, wherein each sub-block corresponds to the Jacobian matrix row of all time channels, and each sub-block is taken as a sub-task;

[0061] The multiple sub-tasks are distributed to corresponding multiple CPU threads;

[0062] Each CPU thread calls a corresponding GPU card to execute the calculation task.

[0063] Further, based on the complete Jacobian matrix, an updated model is obtained, specifically including:

[0064] The gradient and Hessian matrix of the objective function are calculated according to the Jacobian matrix, and the updated model is obtained by solving the normal equation.

[0065] Further, the calculation of the sub-tasks is performed in parallel on the GPUs, specifically including:

[0066] According to the number of meshes divided by the calculation model and the number of observation time channels, the coefficient matrix and multiple right end items of the linear equation system are assembled;

[0067] Forward elimination is performed on the coefficient matrix;

[0068] The backward elimination operation is performed on the multiple right end item vectors, and the unknown solution of the corresponding right end item is calculated;

[0069] The calculated results are assigned to the electric field and magnetic field variables according to the corresponding relationship between the right end items and the unknowns.

[0070] Further, in the process of performing forward elimination on the coefficient matrix:

[0071] The coefficient matrix is pre-decomposed once, the pre-decomposition method is independent of the right end items of the linear system, only the forward elimination processing is performed on the coefficient matrix, the tri-diagonal matrix is simplified into a lower triangular or upper triangular matrix, and the decomposition result is used to solve the equation system corresponding to all right end items.

[0072] Further, the right end items are composed of all observation point data of each detection time channel, and are stored in an interleaved storage format.

[0073] 1. SATEM data acquisition for inversion calculation:

[0074] Semi-Airborne Transient Electromagnetic Method (SATEM) data acquisition adopts ground-air cooperative observation mode. The ground high-power transmitter emits a 1:1 duty cycle trapezoidal wave to generate primary field. The three-component receiving coil system mounted on the unmanned aerial vehicle (UAV) synchronously collects the time-domain decay signal of the vertical magnetic field component rate (d B z / d t ) during the current off period. Wherein B z represents the magnetic field in the z direction. The system uses high-precision GPS synchronous timing technology to ensure the timing synchronization accuracy of the transmitter off and signal acquisition better than 1μs, μs is microsecond. The flight attitude parameters of the UAV are recorded at the same time for motion compensation, and finally the high-quality electromagnetic response data set containing spatial coordinates, time channel information and system parameters is obtained.

[0075] 2. Basic principles of GN inversion method:

[0076] The inversion method deduces the physical parameters of the underground medium (such as resistivity, polarization rate, magnetic permeability, etc.) through observation data. The GN inversion method is a nonlinear optimization algorithm based on second-order approximation, and its core is to reconstruct the underground electrical structure by iteratively solving the minimum value of the objective function. This method converts the inversion problem into a least squares optimization problem, and the objective function contains data fitting terms and model constraint terms. In each iteration, first calculate the forward response of the current model and its Jacobian matrix (sensitivity matrix), then construct the Gauss-Newton equation to solve the model update. Compared with the traditional first-order optimization method, the GN method significantly improves the convergence speed by introducing the approximate calculation of the Hessian matrix. The key factor that determines the computational efficiency of the GN method is the calculation of the sensitivity matrix, which can be expressed as:

[0077] (1)

[0078] wherein d represents the observation data, d i represents the i-th observation data; N r and N m represent the number of observation data and the number of inversion parameters, respectively; represents the Jacobian matrix value of i rows and j columns; m j represents the conductivity of the i-th inversion grid.

[0079] The Jacobian matrix is expressed in the form of a block matrix:

[0080] (2)

[0081] wherein, For the first Jacobian matrix sub-block, For the k-th Jacobian matrix sub-block, For the Nth P A Jacobian matrix sub-block.

[0082] Each sub-block The Jacobian matrix rows corresponding to all time channels are expressed as follows:

[0083] (3)

[0084] This represents the k1-th Jacobian matrix; Denotes the kt-th Jacobian matrix; Indicates the KNth T A Jacobian matrix.

[0085] An operator Q is defined for each sub-block, where Q represents the conversion of the electric field calculated by forward modeling into the electric field at the location of the observed data through curl operator and linear interpolation. B z / d t response.

[0086] The governing equations based on the BEDS-FDTD 3D forward modeling algorithm can be expressed as:

[0087] (4)

[0088] In the formula, vector x represents the solution calculated by forward modeling, and m represents the model parameters. Let m represent the three-dimensional forward modeling operator. Taking the derivative of equation (4) with respect to the model parameter m, we obtain the transpose of the k-th sub-block of the Jacobian matrix:

[0089] (5)

[0090] In the formula, K is the coefficient matrix of the forward equation, and G represents the parameter matrix, specifically:

[0091] (6)

[0092] Where m1 represents the inversion parameters of the first inversion grid, m2 represents the inversion parameters of the second inversion grid, and m Nm This represents the inversion parameters of the Nm-th inversion grid.

[0093] To avoid directly calculating the inverse of matrix K, the adjoint forward modeling method is usually used. Define the adjoint field vector. Then the adjoint forward equation is:

[0094] (7)

[0095] This equation needs to be solved for Nᵣ the right-hand side to obtain the adjoint field vector v. Defining the adjoint field vector v includes the adjoint electric field , the adjoint virtual electric field , and the adjoint magnetic field , the equation (7) is expanded as:

[0096] (8)

[0097] (9)

[0098] (10)

[0099] where η ( η = x , y or z ) denotes a coordinate direction that follows a cyclic shift property (e.g., if η = x , then η +1 and η -1 correspond to y and z , respectively).

[0100] The remaining parameters are defined as:

[0101] , and ; ε denotes permittivity; σ denotes conductivity; Δt denotes time step; μ denotes permeability;

[0102] the superscripts n and n +1 denote the previous time step and the current time step, respectively;

[0103] denotes the first-order difference operator in the η direction; denotes the first-order difference operator in the η-1 direction; denotes the first-order difference operator in the η+1 direction; denotes the second-order difference operator in the η+1 direction; denotes the second-order difference operator in the η-1 direction.

[0104] is an auxiliary field without physical meaning, and denote the adjoint electric field and the adjoint magnetic field, respectively; denotesη direction n+1 the accompanying electric field at the time instant; representing η-1 direction n+1 the accompanying electric field at the time instant; representing η+1 direction n+1 the accompanying magnetic field at the time instant; representing η direction n the accompanying magnetic field at the time instant; representing η direction n+1 the accompanying magnetic field at the time instant; representing η+1 the accompanying virtual electric field of direction; representing η-1 the accompanying virtual electric field of direction; representing η direction n the accompanying electric field at the time instant; representing η-1 direction n the accompanying magnetic field at the time instant; representing the source term at the time instant n+1.

[0105] The Jacobian matrix transpose J T is obtained by formula (7). After the Jacobian matrix is calculated, the model Δm can be updated by solving the normal equation.

[0106] 3. Parallel 3D inversion method based on multiple GPUs

[0107] Figure 1 The flow of the multiple-GPU parallel inversion algorithm of the embodiment is shown.

[0108] First, the necessary initial information for inversion is obtained by reading external files, such as the initial model conductivity distribution, the position of the emission source, the position and time point of the observation data, the grid division rule, etc. Then, the initial model is forward calculated to obtain the root mean square error (RMS), and after calculating the regularization factor required by the current program, the Jacobian matrix calculation task is split into multiple subtasks on the CPU memory, and then multiple subtasks are assigned to multiple CPU threads through the OpenMP command, and then each CPU thread calls the corresponding GPU card to execute the calculation task.

[0109] As Figure 1As shown, J(1) split into is assigned to GPU(1), J(2) split into is assigned to GPU(2), and J(n) split into is assigned to GPU(n), and since there can be at most 8 GPUs installed on each motherboard, there can be at most 8 sub-tasks calculated simultaneously each time. Each GPU opens a kernel computing area, and K1, K2, K n respectively represent the first kernel computing area, the second kernel computing area, and the nth kernel computing area.

[0110] Since different GPUs are independent of each other, after different GPUs perform their respective sub-tasks in parallel, the calculation results are transmitted back to the CPU memory, and other sub-tasks are executed for calculation until all tasks are completed, and the results of all sub-tasks are spliced into a complete Jacobian matrix.

[0111] According to the complete Jacobian matrix, the gradient and Hessian matrix of the objective function are calculated, and the solution of the equation is obtained. Then the updated model is obtained. After that, forward calculation is performed again to obtain RMS and judge whether the preset condition is met. If not, the regularization factor is modified and the Jacobian matrix is calculated repeatedly until the RMS meets the calculation threshold. Finally, the current conductivity model is output, which is the three-dimensional inversion result.

[0112] It can be seen that solving the Jacobian matrix is the core part of the algorithm, in order to improve the calculation efficiency of the Jacobian matrix, a multi-GPU parallel strategy is adopted:

[0113] OpenMP technology is used to control CPU multi-thread management of multi-GPU devices to realize parallel calculation of block matrices (equation (2)); OpenACC and CUDA parallel calculation of sub-matrices (equation (3)) are combined on a single GPU device; and finally the assembly of the Jacobian matrix J is completed on the CPU side. OpenACC and CUDA are operation command sets that can be specified on GPU devices.

[0114] According to equation (7), multiple right-hand side (RHS) tridiagonal linear systems are required to be solved for each iteration of the forward. Therefore, a multi-GPU parallel computing scheme is proposed to speed up the solving process. When considering the time channel data of all observation points, a tridiagonal linear equation group containing multiple RHS as shown in equation (8) needs to be solved. Figure 2 In equation (8), the tridiagonal matrix on the left side is the coefficient matrix, and the black dots represent the positions with numerical values; the matrix on the right side represents the multiple right-hand sides, and each column uses different depth of gray dots to represent the numerical values corresponding to the pseudo-source terms; and the matrix in the middle represents the unknown terms to be solved. Figure 2

[0115] ​The traditional cuThomasBatch solver cannot solve linear equations with multiple right-hand side terms, so it can only copy multiple coefficient matrices according to the number of RHS and assemble multiple linear equations into a large linear system for solving, which will seriously reduce the calculation efficiency.

[0116] 4. A method for solving multiple-RHS linear equations based on the Thomas Batch algorithm

[0117] To solve the above problems, a method for solving multiple right-hand side tridiagonal linear equations based on the Thomas algorithm is applied. The Thomas algorithm is essentially a simplified form of the Gaussian elimination method, which includes a forward pass and a backward pass. The flowchart of the Thomas algorithm for solving tridiagonal matrices is shown in FIG. 1. Figure 3 In the forward pass stage, the tridiagonal equation system is transformed into an upper triangular form, and then the upper triangular non-diagonal elements are eliminated by the backward pass, and finally the diagonalized equation system is obtained, so that the unknowns can be directly solved.

[0118] Figure 3 The flowchart of the Thomas algorithm for solving tridiagonal matrices is shown in FIG. 1. Specifically, Figure 3 (a) in FIG. 1 shows a tridiagonal linear system, Figure 3 (b) in FIG. 1 shows the forward pass process, Figure 3 (c) in FIG. 1 shows the backward pass process, where the black circles in the band matrix represent non-zero terms. The forward pass eliminates the lower diagonal elements, and the backward pass eliminates the upper diagonal elements. For simplicity, the RHS and unknown vectors are omitted.

[0119] The specific implementation process of the Thomas algorithm is shown in Algorithm 1 as follows. First, set the variables a, b, c, and d to store the sparse matrix tridiagonal elements and the right-hand side terms, where the subscript s represents the s-th loop, the subscript s-1 represents the s-1-th loop, and the subscript s+1 represents the s+1-th loop; then calculate u1=d1 / b1 and c1=c1 / b1 (lines 1-2). u1 represents an intermediate variable used to store data in the first loop, d1 represents the right-hand side in the first loop, a1 represents the upper triangular element in the first loop, b1 represents the diagonal element in the first loop, and c1 represents the lower triangular element in the first loop.

[0120] By forward elimination, calculate ω ← 1 / ( b s ← a s c s - 1 ), u s ←ω ( d s ← a s u s - 1 ), c s ← ω c s , where s is the cycle indicator (lines 3-7). ω This represents an intermediate variable used to store data. b s This represents the diagonal element in the s-th iteration of the loop. a s This represents the upper triangular element in the s-th cycle. c s – 1 This represents the lower triangular element in the (s-1)th iteration of the cycle. d s This represents the right-hand term in the s-th iteration of the cycle. u s – 1 This represents the intermediate variable used to store data in the (s-1)th iteration of the loop.

[0121] Then perform backward elimination. ω ← 1 / ( b s ← a s c s - 1 )and u s ← u s - c s u s + 1 (Lines 8-11), and obtain the unknown to be solved. u s This represents the intermediate variable used to store data in the s-th iteration of the loop. u s + 1 This represents the intermediate variable used to store data in the (s+1)th iteration of the loop.

[0122]

[0123] For linear equation systems with multiple right-hand sides (RHS), although the coefficient matrix decomposition requires adjustments for each right-hand side term, the resulting matrix can be shared by multiple right-hand side terms. The algorithm implementation steps are as follows:

[0124] (1) Initialize variables: Initialize the variable space of electric field, magnetic field, coefficient matrix, right-hand side, etc.; initialize the handle space of CUDA functions and initialize the variable space required in all GPU memory.

[0125] (2) Open up the CPU multi-threaded computing area: decompose the Jacobian matrix into multiple sub-tasks, and distribute them to multi-threaded execution calculation according to the number of observation points through the OpenMP command.

[0126] (3) Open all GPU computing areas: accept CPU multi-threaded control of all GPU execution sub-tasks through the CUDA command.

[0127] (4) Each GPU opens kernel computing area 1: according to the number of grids and the number of observation time channels, assemble the coefficient matrix and multiple RHS of the linear equation group.

[0128] (5) Each GPU opens kernel computing area 2: performs a pre-decomposition operation on the coefficient matrix, that is, performs forward elimination.

[0129] (6) Each GPU opens kernel computing area 3: performs backward elimination operation on the multiple RHS vector, and calculates the unknown solution corresponding to the RHS.

[0130] (7) Each GPU opens kernel computing area 4: assigns the calculated results to the electric field and magnetic field variables.

[0131] (8) Release the GPU computing area: copy the electric field and magnetic field variables from the video memory back to the CPU memory.

[0132] (9) Release the CPU's non-main thread: assemble the calculation results of all sub-thread sub-tasks into the Jacobian matrix.

[0133] The embodiment of the application proposes a multi-GPU parallel sensitivity matrix solving scheme based on the BEDS-FDTD forward framework. According to the characteristics of the GN optimization method, the task scheduling and load balancing of multiple CPU cores are realized through OpenMP, and efficient parallel computing on the GPU end is realized by combining OpenACC technology, and a heterogeneous parallel computing architecture of CPU master control and multi-GPU cooperation is constructed.

[0134] In specific implementation, the CPU is responsible for overall task decomposition, memory management and device synchronization, while multiple GPUs parallelly process the sub-matrix calculation tasks allocated respectively. This scheme fully utilizes the parallel computing capability of modern computing devices, provides an efficient and scalable computing framework for time-domain electromagnetic inversion, significantly improves the calculation efficiency of the sensitivity matrix in the GN inversion, and provides a new technical approach for solving large-scale electromagnetic inversion problems.

[0135] For the adjoint forward control equation, the embodiment of the application innovatively puts forward an improved ThomasBatch pre-decomposition batch solution method. According to the characteristics of the adjoint forward equation based on the BEDS-FDTD forward algorithm, the coefficient matrix corresponding to multiple right-hand side terms (RHS) has the same structural characteristics, and the traditional solving process is improved to a batch processing mode. Specifically, the method first performs one-time pre-decomposition on the coefficient matrix, and then uses the decomposition result to solve all the equation groups corresponding to the right-hand side terms at the same time, avoiding the repeated matrix decomposition operation in the traditional cuThomasBatch solver. Especially for the characteristics of dense observation data in the semi-airborne transient electromagnetic (SATEM) method, the method further optimizes the memory access mode and the calculation process, significantly improving the solving efficiency. This innovation not only provides an efficient solution scheme for SATEM three-dimensional inversion, but also can be applied to other geophysical problems that need to solve multiple right-hand side linear systems, and has important theoretical value and practical application significance.

[0136] To verify the performance of the algorithm proposed in this embodiment, the sensitivity matrix J is calculated and the relevant calculation time is extracted (the initial time cost such as CPU / GPU memory allocation, pseudo-source setting, etc. is ignored). The experimental platform is configured as: 8 NVIDIA Tesla A100 (80GB memory), Intel(R) Xeon(R) Platinum 8350C processor (128 cores) and 1003GB memory. By defining the speedup ratio (formula 11), the algorithm efficiency is compared.

[0137] (11)

[0138] In the formula, R represents the calculation speed, and the unit is s. The speedup ratio is represented by The calculation time of the ThomasBatch algorithm is represented by The calculation time of the cuThomasBatch algorithm is represented by

[0139] Figure 4 The speedup effect of the ThomasBatch algorithm relative to the cuThomasBatch algorithm is shown. Figure 4 The horizontal coordinate in the figure represents the number of right-hand side terms, and the vertical coordinate represents the size of the linear equation group. Different speedup ratios are represented by different degrees of blue, where the deeper the color degree, the greater the speedup ratio. The experimental results show that the speedup ratio of the ThomasBatch algorithm is always higher than 2.71, proving that its calculation efficiency is significantly better than that of the cuThomasBatch algorithm based on the sensitivity matrix.

[0140] From Figure 4It can be seen that all the algorithms achieve obvious acceleration, but their performance characteristics are different in different cases. The acceleration effect of the ThomasBatch algorithm gradually decreases with the increase of the number of grids.

[0141] Embodiment two

[0142] The embodiment discloses a multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion system.

[0143] The multi-GPU parallel-based semi-airborne transient electromagnetic fast inversion system comprises:

[0144] A data acquisition module is configured to acquire necessary initial information for semi-airborne transient electromagnetic inversion.

[0145] A model establishment module is configured to establish an initial model and determine a regularization factor of the initial model.

[0146] A Jacobian matrix determination module is configured to determine a calculation task of a Jacobian matrix based on the necessary initial information and the initial model.

[0147] A task allocation module is configured to split the calculation task of the Jacobian matrix into multiple subtasks on a CPU memory, transmit the multiple subtasks to corresponding GPUs respectively, and perform parallel calculation of the subtasks on the GPUs.

[0148] A parallel execution module is configured to, after the GPUs perform the respective subtasks in parallel, transmit the calculation results back to the CPU memory, and perform calculation of other subtasks until the execution of all the subtasks is completed.

[0149] A return splicing module is configured to splice the results of all the subtasks into a complete Jacobian matrix on the CPU memory.

[0150] A model updating module is configured to obtain an updated model based on the complete Jacobian matrix.

[0151] A judgment loop module is configured to perform forward calculation on the updated model to obtain a root mean square error, judge whether the root mean square error meets a threshold value, modify the regularization factor and repeatedly calculate the Jacobian matrix if the root mean square error does not meet the threshold value, and obtain a final model as a three-dimensional inversion result until the root mean square error meets the threshold value.

[0152] Those skilled in the art should understand that the modules or steps of the present application described above can be realized by general computer devices, or alternatively, they can be realized by program codes executable by the computer devices, so that they can be stored in the storage devices and executed by the computer devices, or they can be respectively manufactured into individual integrated circuit modules, or a plurality of modules or steps among them can be manufactured into a single integrated circuit module. The present application is not limited to any specific combination of hardware and software.

[0153] The specific embodiments of the present application described above with reference to the drawings are not intended to limit the scope of protection of the present application. Those skilled in the art should understand that various modifications or variations made on the basis of the technical solutions of the present application without creative labor are still within the scope of protection of the present application.

Claims

1. A fast inversion method for semi-airborne transient electromagnetic data based on multi-GPU parallel processing, characterized in that: Includes the following steps: To obtain the necessary initial information for semi-airborne transient electromagnetic inversion; Establish an initial model and determine its regularization factor; Based on the necessary initial information and initial model, the task of calculating the Jacobian matrix is ​​determined; The Jacobian matrix calculation task is split into multiple subtasks on the CPU memory, and these subtasks are transferred to their respective GPUs for parallel computation. Specifically, this includes: The Jacobian matrix is ​​expressed in block matrix form, where each sub-block corresponds to a row of the Jacobian matrix for all time channels, and each sub-block is treated as a subtask. Assign multiple subtasks to corresponding CPU threads; Each CPU thread calls the corresponding GPU to perform computing tasks; Based on the number of grids and observation time channels in the computational model, assemble the coefficient matrix and multiple right-hand terms of the linear equation system; Perform forward elimination on the coefficient matrix; perform backward elimination on the vector of multiple right-hand terms and calculate the unknown solutions for the corresponding right-hand terms; assign the calculation results to the electric and magnetic field variables according to the correspondence between the right-hand terms and the unknowns; perform a one-time pre-decomposition of the right-hand terms of the coefficient matrix that is decoupled from the linear system, and perform forward elimination only on the coefficient matrix to simplify the tridiagonal matrix into a lower triangular or upper triangular matrix to obtain the decomposition results; use the decomposition results to solve the system of equations corresponding to all right-hand terms simultaneously; The coefficient matrices corresponding to multiple right-hand terms have the same structural characteristics, which improves the traditional process of solving them one by one into a batch processing mode. Specifically, the coefficient matrix is ​​first pre-decomposed, and then the decomposition results are used to solve the system of equations corresponding to all right-hand terms at the same time. After the GPU completes its respective subtask in parallel, it sends the calculation results back to the CPU memory and then takes over other subtasks to perform calculations until all subtasks have been completed. Concatenate the results of all subtasks into a complete Jacobian matrix in the order of their splitting on the CPU memory. The updated model is obtained based on the complete Jacobian matrix; The updated model is subjected to forward modeling to obtain the root mean square error. It is then determined whether the root mean square error meets the threshold. If it does not, the regularization factor is modified and the Jacobian matrix is ​​repeatedly calculated in a loop until the root mean square error meets the threshold. The final model obtained is the three-dimensional inversion result.

2. The fast inversion method for semi-airborne transient electromagnetic events based on multi-GPU parallel processing as described in claim 1, characterized in that, The necessary initial information for the semi-airborne transient electromagnetic inversion includes the conductivity distribution of the initial model, the location and size of the emission source, the location and time of the observation data, the grid size, the boundary conditions, and the grid partitioning rules.

3. The fast inversion method for semi-airborne transient electromagnetic events based on multi-GPU parallel processing as described in claim 1, characterized in that, The initial model represents a meshed model containing a pre-defined distribution of electrical parameters before calculation; the updated model represents a meshed model containing a distribution of electrical parameters calculated during the inversion process; and the final model represents a meshed model containing the final obtained distribution of electrical parameters.

4. The fast inversion method for semi-airborne transient electromagnetic events based on multi-GPU parallel processing as described in claim 1, characterized in that, Based on the necessary initial information and initial model, the task of calculating the Jacobian matrix is ​​determined as follows: The Jacobian matrix is ​​represented as: ； In the formula, d Represents observation data, N r and N m These represent the number of observation data and the number of inversion parameters, respectively. Let represent the Jacobian matrix in row i and column j; d i This represents the i-th observation data; m j Let represent the conductivity of the i-th inversion grid.

5. The fast inversion method for semi-airborne transient electromagnetic events based on multi-GPU parallel processing as described in claim 1, characterized in that, Based on the complete Jacobian matrix, the updated model is obtained, specifically including: The gradient of the objective function and the Hessian matrix are calculated using the Jacobian matrix, and the updated model can be obtained by solving the normal equations.

6. The fast inversion method for semi-airborne transient electromagnetic events based on multi-GPU parallel processing as described in claim 1, characterized in that, The right-hand item consists of all observation point data from each detection time channel, and is stored in an interleaved storage format.

7. A semi-airborne transient electromagnetic fast inversion system based on multi-GPU parallel processing, characterized in that, include: The data acquisition module is configured to acquire the necessary initial information for semi-airborne transient electromagnetic inversion. The model building module is configured to: build an initial model and determine the regularization factor of the initial model; The Jacobian matrix determination module is configured to: determine the task of calculating the Jacobian matrix based on necessary initial information and an initial model; The task allocation module is configured to: split the Jacobian matrix calculation task into multiple subtasks on the CPU memory, transfer the multiple subtasks to the corresponding GPUs respectively, and execute the calculation of the subtasks in parallel on the GPUs; The parallel execution module is configured to: after the GPU completes its respective subtask in parallel, transfer the calculation results back to the CPU memory, and then retrieve other subtasks for calculation until all subtasks have been completed, specifically including: The Jacobian matrix is ​​expressed in block matrix form, where each sub-block corresponds to a row of the Jacobian matrix for all time channels, and each sub-block is treated as a subtask. Assign multiple subtasks to corresponding CPU threads; Each CPU thread calls the corresponding GPU to perform computing tasks; Based on the number of grids and observation time channels in the computational model, assemble the coefficient matrix and multiple right-hand terms of the linear equation system; Perform forward elimination on the coefficient matrix; perform backward elimination on the vector of multiple right-hand terms and calculate the unknown solutions for the corresponding right-hand terms; assign the calculation results to the electric and magnetic field variables according to the correspondence between the right-hand terms and the unknowns; perform a one-time pre-decomposition of the right-hand terms of the coefficient matrix that is decoupled from the linear system, and perform forward elimination only on the coefficient matrix to simplify the tridiagonal matrix into a lower triangular or upper triangular matrix to obtain the decomposition results; use the decomposition results to solve the system of equations corresponding to all right-hand terms simultaneously; The coefficient matrices corresponding to multiple right-hand terms have the same structural characteristics, which improves the traditional process of solving them one by one into a batch processing mode. Specifically, the coefficient matrix is ​​first pre-decomposed, and then the decomposition results are used to solve the system of equations corresponding to all right-hand terms at the same time. The returned concatenation module is configured to concatenate the results of all subtasks into a complete Jacobian matrix on the CPU memory in the order of splitting. The model update module is configured to obtain the updated model based on the complete Jacobian matrix; The judgment loop module is configured to: perform forward calculation on the updated model to obtain the root mean square error, determine whether the root mean square error meets the threshold, if not, modify the regularization factor and repeat the loop calculation of the Jacobian matrix until the root mean square error meets the threshold, and the final model obtained is the three-dimensional inversion result.

Images