Microseismic reverse time imaging method, system and device based on BF16 mixing precision and medium

Abstract

The invention belongs to the technical field of oil and gas geophysical exploration, and relates to a microseism reverse time imaging method, system and device based on BF16 mixing precision and a medium. According to the method, a BF16 and FP32 mixed storage format is used, a heterogeneous parallel computing framework composed of a CPU and a GPU is adopted to execute a microseism reverse time imaging method based on BF16 mixed precision, a BF16 mixed precision computing strategy and CPU + GPU heterogeneous parallelism are introduced, the computing speed of seismic source reverse time imaging can be increased by several times, and the computing speed of seismic source reverse time imaging can be increased by several times. A BF16 format is adopted and assisted by a high-precision accumulator and an FP32 main copy storage strategy, so that the performance advantage brought by low-precision calculation is enjoyed, the numerical stability under long-term iteration and the physical reliability of a final result are ensured, and efficient calculation and low-storage use of micro-seismic reverse time imaging are realized.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas geophysical exploration technology, and relates to a microseismic reverse time imaging method, system, equipment and medium based on BF16 mixed accuracy. Background Technology

[0002] Hydraulic fracturing is a core technology for the commercial exploitation of shale gas. It involves injecting high-pressure fluid into the ground to create a network of fractures, thereby increasing reservoir permeability. During this process, the fracturing of the rock formations induces numerous weak seismic events, known as "microseismic events." Accurate monitoring of the spatiotemporal distribution of these microseismic events is crucial for assessing fracturing effectiveness, depicting fracture network morphology, and predicting production capacity.

[0003] Time-Reverse Imaging (TRI) is a cutting-edge technology in the field of microseismic monitoring. Based on wave equation theory, this technology reverses the temporal order of complete microseismic waveform data recorded by ground stations, using it as a "wave source" that propagates back into the subsurface medium. By applying specific imaging conditions, the energy of different types of seismic waves (such as P-waves and S-waves) is focused at the actual source location and origin time, thus achieving precise localization of microseismic events. Compared to traditional localization methods that rely on time-based image acquisition, TRI eliminates the need for manual intervention in image acquisition, is more adaptable to low signal-to-noise ratio data, and theoretically can handle arbitrarily complex subsurface media.

[0004] However, existing seismic source reverse-time imaging technology still faces significant challenges in practical applications, primarily due to the bottleneck of computational resources. Seismic source reverse-time imaging is a computationally and storage-intensive process. It involves finely meshing three-dimensional space and iteratively solving the wave equation over thousands or even tens of thousands of time steps, requiring enormous computational resources and memory (especially GPU memory). This bottleneck not only makes its computational cost high, its cycle lengthy, and its computational efficiency low, greatly limiting its widespread application in production, but also makes it difficult to efficiently integrate with artificial intelligence technologies such as deep learning, which require massive amounts of data for model training, thus failing to meet the demands of deep learning and other AI technologies for rapid processing of massive amounts of data. Summary of the Invention

[0005] The purpose of this invention is to provide a microseismic reverse-time imaging method, system, device, and medium based on BF16 mixed precision, in order to solve the technical problems of low computational efficiency and inability to meet the needs of deep learning and other artificial intelligence technologies for rapid processing of massive amounts of data in practical applications of existing source reverse-time imaging technology.

[0006] To achieve the above objectives, the present invention employs the following technical solution: In a first aspect, the present invention provides a microseismic reverse-time imaging method based on BF16 hybrid precision, using a hybrid storage format of BF16 and FP32, and executing it using a heterogeneous parallel computing framework composed of CPU and GPU, including the following steps: A three-dimensional arbitrary curve physical mesh is generated based on surface elevation data, and the three-dimensional arbitrary curve physical mesh is stored in FP32 format. The microseismic data recorded by the station was read in FP32 format, and the microseismic data was reversed on the time axis to obtain the reverse time reconstructed wavefield; Using the reconstructed wavefield in reverse time as the boundary source term, on a three-dimensional arbitrary curved physical mesh, the finite difference method of the curved mesh in BF16 format is used to perform time iteration core calculations on all time steps on the GPU through CUDA, and the first-order velocity-stress equations are solved iteratively to obtain the wavefield master copy. The P-wave and S-wave acceleration fields are decoupled and separated from the primary and secondary copies of the wave field using the BF16 scheme to obtain the P-wave acceleration field and the S-wave acceleration field. Based on the longitudinal wave acceleration field and the transverse wave acceleration field, and based on the imaging condition formula, the imaging values ​​at the grid points are obtained in BF16 format. The imaging values ​​at the grid points of all time steps are aggregated, and the aggregation result is transmitted back to the CPU in FP32 format to construct a four-dimensional imaging data volume. The maximum value in the four-dimensional imaging data volume is searched globally. The maximum value is the location of the microseismic event, and the time corresponding to the maximum value is the time when the microseismic event occurred.

[0007] Secondly, the present invention provides a microseismic reverse-time imaging system based on BF16 hybrid precision, comprising: A three-dimensional arbitrary curve physical mesh acquisition module is used to generate a three-dimensional arbitrary curve physical mesh based on surface elevation data and store the three-dimensional arbitrary curve physical mesh in FP32 format. The reverse-time reconstructed wavefield acquisition module is used to read the microseismic data recorded by the station in FP32 format, reverse the microseismic data on the time axis, and obtain the reverse-time reconstructed wavefield. The elastic wave propagation result acquisition module is used to reconstruct the wave field in reverse time as the boundary source term. On a three-dimensional arbitrary curve physical grid, the time iterative core calculation is performed on the GPU using CUDA in the BF16 format and the curve grid finite difference method to perform time iterative core calculations on all time steps. The first-order velocity-stress equations are solved iteratively to obtain the wave field master copy. The acceleration field decomposition module is used to decouple and separate the longitudinal and transverse waves in the main copy of the wave field in BF16 format to obtain the longitudinal wave acceleration field and the transverse wave acceleration field. The imaging value acquisition module at the grid point is used to acquire the BF16 format imaging value at the grid point based on the longitudinal wave acceleration field and the transverse wave acceleration field and the imaging condition formula. The four-dimensional imaging data volume acquisition module is used to aggregate the imaging values ​​at the grid points at all time steps, and transmit the aggregation result back to the CPU in FP32 format to construct the four-dimensional imaging data volume. The microseismic event confirmation module is used to globally search for the maximum value in the four-dimensional imaging data volume. The maximum value is the location of the microseismic event, and the time corresponding to the maximum value is the time when the microseismic event occurred.

[0008] Thirdly, the present invention provides an electronic device, comprising: a processor; a memory for storing computer program instructions; and steps for implementing a microseismic reverse-time imaging method based on BF16 mixed precision when executing the computer program.

[0009] Fourthly, the present invention provides a storage medium storing computer program instructions, which are loaded and executed by a processor, wherein the processor performs a microseismic reverse-time imaging method based on BF16 mixed precision.

[0010] Compared with the prior art, the present invention has the following beneficial effects: This invention generates a three-dimensional arbitrary curve physical grid based on surface elevation data, which can more accurately fit the complex terrain. The three-dimensional arbitrary curve physical grid is stored in FP32 format, providing a foundation for subsequent fine simulation of microseismic wave propagation and helping to ensure the accuracy of geological simulation. Microseismic data recorded by stations is read in FP32 format, and the data is reversed on the time axis. Through time reversal, the actually recorded microseismic waveform data is transformed into "wave source" data that can be used for reverse-time imaging. This allows the wavefield to propagate back to the subsurface medium, laying the foundation for energy focusing at the actual source location and origin time, thereby achieving precise localization of microseismic events. Using the finite difference method of the curve grid, time iteration core calculations are performed on all time steps on the GPU using CUDA, iteratively solving the first-order velocity-stress equations to obtain a master copy of the wavefield. By solving the first-order velocity-stress equations in reverse time, the propagation state of elastic waves at different times and locations can be obtained, providing rich wavefield data for further analysis of microseismic events. Using the BF16 format and the finite difference method with curved meshes to adapt to complex meshes, parallel computation is performed on the GPU via CUDA. Combined with BF16 mixed precision, this significantly improves computational speed, reduces memory and GPU memory usage, and resolves computational resource bottlenecks. The BF16 format is used to decouple P-waves and S-waves from the primary and secondary wave copies of the wavefield, obtaining P-wave acceleration fields and S-wave acceleration fields. This helps to more accurately understand the generation mechanism and propagation process of microseismic events, providing more precise wavefield information for subsequent acquisition of imaging values ​​based on imaging condition formulas. Based on the P-wave and S-wave acceleration fields and the imaging condition formulas, BF16 format imaging values ​​at grid points can be obtained, allowing for the preliminary determination of potential microseismic event areas in three-dimensional space, providing crucial evidence for the final precise location of microseismic events. The imaging values ​​at the grid points across all time steps are aggregated, and the aggregation results are transmitted back to the CPU in FP32 format to construct a four-dimensional imaging data volume, providing rich data support for subsequent in-depth analysis and assessment of microseismic events. The invention achieves precise localization of microseismic events by globally searching for the maximum value in the four-dimensional imaging data volume, where the maximum value represents the location of the microseismic event and the corresponding time represents the occurrence time of the microseismic event. By introducing a BF16 mixed-precision computing strategy and CPU+GPU heterogeneous parallelism, the computation speed of source reverse-time imaging can be increased several times. Employing the BF16 format supplemented by a high-precision accumulator and an FP32 master-replica storage strategy, it enjoys the performance advantages of low-precision computing while ensuring numerical stability and physical reliability of the final results over long-term iterations, thus achieving efficient computation and low storage usage for microseismic reverse-time imaging.

[0011] The system of this invention includes a three-dimensional arbitrary curve physical mesh acquisition module, a reverse-time reconstruction wavefield acquisition module, an elastic wave propagation result acquisition module, an acceleration field decomposition module, an imaging value acquisition module at grid points, a four-dimensional imaging data volume acquisition module, and a microseismic event confirmation module. The three-dimensional arbitrary curve physical mesh acquisition module generates a three-dimensional arbitrary curve physical mesh based on surface elevation data and stores the three-dimensional arbitrary curve physical mesh in FP32 format. The reverse-time reconstruction wavefield acquisition module reads microseismic data recorded by stations in FP32 format, reverses the microseismic data on the time axis, and obtains a reverse-time reconstructed wavefield. The elastic wave propagation result acquisition module uses the reverse-time reconstructed wavefield as a boundary source term and, on the three-dimensional arbitrary curve physical mesh, performs time iteration core calculations on all time steps using the curve mesh finite difference method in BF16 format on the GPU via CUDA, iteratively solving the first-order velocity-stress equations to obtain a master copy of the wavefield. The acceleration field decomposition module... The decoupling module is used to decouple and separate the P-wave and S-wave in the wavefield master copy in BF16 format to obtain the P-wave acceleration field and the S-wave acceleration field. The grid point imaging value acquisition module is used to obtain the BF16 format imaging value at the grid point based on the P-wave acceleration field and the S-wave acceleration field and the imaging condition formula. The four-dimensional imaging data volume acquisition module is used to aggregate the imaging values ​​at the grid point at all time steps and transmit the aggregation result back to the CPU in FP32 format to construct the four-dimensional imaging data volume. The microseismic event confirmation module is used to globally search for the maximum value in the four-dimensional imaging data volume. The maximum value is the location of the microseismic event, and the time corresponding to the maximum value is the time of occurrence of the microseismic event. The various modules work together to achieve efficient computation and low storage usage for microseismic reverse-time imaging.

[0012] The electronic device and storage medium of this invention can also achieve efficient computation and low storage usage for microseismic reverse time imaging. Attached Figure Description

[0013] Figure 1 This is a flowchart of a method according to an embodiment of the present invention; Figure 2 This is a system block diagram of an embodiment of the present invention; Figure 3 This is a flowchart of the CPU+GPU heterogeneous parallel acceleration algorithm for seismic source reverse-time imaging based on BF16 hybrid precision, according to an embodiment of the present invention. Figure 4 A schematic diagram illustrating the storage formats of floating-point numbers of different precisions (FP32, FP16, BF16) in a computer; Figure 5 This is a P-wave velocity model diagram of a three-dimensional complex medium model according to an embodiment of the present invention; Figure 6The results of microseismic reverse-time imaging localization tests on a three-dimensional complex medium model according to an embodiment of the present invention are shown below. Figure 6 a represents the source location results based on BF16 hybrid precision reverse-time imaging. Figure 6 b represents the source location result from reverse-time imaging with standard FP32 precision; Figure 7 The image is an error histogram comparing the BF16 mixed-precision imaging results and the conventional FP32 precision imaging results of a three-dimensional complex medium model according to an embodiment of the present invention. Figure 8 This is a comparison chart of computational efficiency and resource consumption between the BF16 mixed-precision method and the conventional FP32 precision method for three-dimensional complex media models according to embodiments of the present invention. Figure 8 'a' represents the comparison of calculation times. Figure 8 b represents a comparison of peak video memory usage; Figure 9 This is a velocity model and topographic relief data for a hydraulic fracturing operation area in a shale gas extraction project. Figure 9 a represents the P-wave velocity model. Figure 9 b represents the S-wave velocity model; Figure 10 This is microseismic monitoring data for hydraulic fracturing in shale gas extraction in a certain region, as described in this embodiment of the invention. Figure 11 This is the source reversal time imaging localization result of microseismic monitoring data from hydraulic fracturing for shale gas extraction in a certain region, according to an embodiment of the present invention. Figure 11 a represents the source location results based on BF16 hybrid precision reverse-time imaging. Figure 11 b represents the source location result from reverse-time imaging with standard FP32 precision; Figure 12 The image is an error histogram of the BF16 mixed precision imaging results and the conventional FP32 precision imaging results of the measured microseismic data in this embodiment of the invention. Figure 13 This chart compares the computational efficiency and resource consumption of the BF16 hybrid accuracy method and the conventional FP32 accuracy method for measured microseismic data according to an embodiment of the present invention. Figure 13 'a' represents the comparison of calculation times. Figure 13 b represents the peak video memory usage comparison. Detailed Implementation

[0014] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0015] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0016] The present invention will now be described in further detail with reference to the accompanying drawings: Example 1: See Figure 1 This embodiment discloses a microseismic reverse-time imaging method based on BF16 mixed precision. It uses a mixed storage format of BF16 (16-bit floating-point) and FP32 (32-bit single-precision floating-point), and is executed using a heterogeneous parallel computing framework composed of a CPU (Central Processing Unit) and a GPU (Graphics Processing Unit). The method includes the following steps: S1 generates a three-dimensional arbitrary curve physical grid based on surface elevation data, which can more accurately fit the actual complex terrain and make the subsequent microseismic monitoring simulation more in line with the real surface environment. The three-dimensional arbitrary curve physical grid is stored in FP32 format. The FP32 format storage ensures the accuracy of the grid data and provides a foundation for the subsequent fine simulation of microseismic wave propagation, overcoming the shortcomings of traditional methods in geological simulation accuracy.

[0017] In one embodiment of the present invention, a three-dimensional arbitrary curve physical grid is generated based on surface elevation data, and a mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system is established based on the three-dimensional arbitrary curve physical grid.

[0018] S2 reads the microseismic data recorded by the station in FP32 format, reverses the microseismic data on the time axis, and obtains the reverse-time reconstructed wavefield. Through the time reversal operation, the actually recorded microseismic waveform data is transformed into "wave source" data that can be used for reverse-time imaging, enabling the wavefield to be propagated back to the subsurface medium in reverse. This lays the foundation for energy focusing at the true source location and origin time, thereby achieving precise localization of microseismic events.

[0019] In this embodiment of the invention, the formula for obtaining the reverse-time reconstructed wavefield is as follows:

[0020] in, From the position of the station to any point in space Green's function; This is the total duration of the data; Represents convolution; To reconstruct the wave field in reverse time; Spatial location coordinates; It is a time variable; Location of the station; Component data recorded by the station.

[0021] S3 uses a reverse-time reconstructed wavefield as the boundary source term. On a 3D arbitrary curved physical mesh, it employs the curved mesh finite difference method in BF16 format, performing time-iterative core calculations on the GPU using CUDA (Compute Unified Device Architecture) for all time steps, iteratively solving the first-order velocity-stress equations to obtain a master copy of the wavefield. The curved mesh finite difference method has high accuracy and can accurately simulate the propagation process of elastic waves on a 3D arbitrary curved physical mesh. By solving the first-order velocity-stress equations in reverse time, the propagation state of elastic waves at different times and locations can be obtained, providing rich wavefield data for further analysis of microseismic events. The use of the curved mesh finite difference method adapts to complex meshes, and the parallel computation on the GPU using CUDA, combined with BF16 mixed precision, significantly improves computation speed, reduces memory and GPU memory usage, and solves the computational resource bottleneck problem.

[0022] In this embodiment of the invention, the step of performing time-iteration core calculations for all time steps, iteratively solving the first-order velocity-stress equations to obtain the wavefield master copy includes: A mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system is established based on the aforementioned three-dimensional arbitrary curve physical mesh; By using the mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system, the first-order velocity-stress equations are transformed into the curve grid coordinate system to obtain the transformed first-order velocity-stress equations. Based on the transformed first-order velocity-stress equations, at each time step, the wave field value is read from the FP32 global memory, the wave field value is converted to BF16 format in the FP32 local accumulator, multiplication is performed and the value is accumulated to obtain the accumulator result, and then the accumulator result is converted back to FP32 format and written back to the global memory as the main copy of the wave field.

[0023] The first-order velocity-stress equations are as follows:

[0024] in, Density; For speed Quantity; For speed Quantity; For speed Quantity; These are the components of the stress tensor; and For Lamé parameters; for Coordinates of direction; for Coordinates of direction; for Coordinates of direction.

[0025] S4. Using the BF16 format, the primary and secondary wave fields of the wavefield are decoupled and separated into P-wave and S-wave acceleration fields to obtain the P-wave acceleration field and the S-wave acceleration field. P-waves and S-waves have different propagation characteristics and response modes to the subsurface medium. Decoupling and separating the total acceleration wavefield into P-wave acceleration fields and S-wave acceleration fields allows for the separate analysis of the roles and propagation laws of P-waves and S-waves in microseismic events. This helps to more accurately understand the generation mechanism and propagation process of microseismic events, and provides more precise wavefield information for subsequent acquisition of imaging values ​​based on imaging condition formulas.

[0026] S5. Based on the P-wave and S-wave acceleration fields and the imaging condition formula, obtain the BF16 format imaging values ​​at the grid points. The imaging condition formula can convert the energy information of P-waves and S-waves into imaging values ​​at the grid points. These imaging values ​​reflect the probability of a seismic source at different locations. By calculating the imaging values, the area where microseismic events may occur can be preliminarily determined in three-dimensional space, providing a key basis for the final precise location of microseismic events. Specifically: The energy propagation direction vector of the longitudinal wave acceleration field is obtained from the longitudinal wave acceleration field. The formula for obtaining the energy propagation direction vector of the longitudinal wave acceleration field is:

[0027] Obtain the energy propagation direction vector of the shear wave acceleration field from the shear wave acceleration field; The formula for obtaining the energy propagation direction vector of the transverse wave acceleration field is:

[0028] in, For the permutation tensor; This is the energy propagation direction vector of the longitudinal wave acceleration field; This is the energy propagation direction vector of the transverse wave acceleration field; Let be the divergence of the velocity field; For longitudinal wave acceleration wave field; and For the permutation tensor; For transverse wave acceleration wave field; Spatial derivative of the velocity field.

[0029] Based on the energy propagation direction vectors of the longitudinal wave acceleration field and the transverse wave acceleration field, and using the imaging condition formula, the imaging values ​​at the grid points are obtained. The imaging condition formula is as follows:

[0030] in, For grid points At any moment The imaging value; This is the energy propagation direction vector of the longitudinal wave acceleration field; This is the energy propagation direction vector of the transverse wave acceleration field.

[0031] S6. Aggregate the imaging values ​​at the grid points of all time steps, and transmit the aggregation result back to the CPU in FP32 format to construct a four-dimensional imaging data volume. The multi-dimensional data volume can more comprehensively and accurately describe the spatiotemporal characteristics of microseismic events, providing rich data support for subsequent in-depth analysis and assessment of microseismic events.

[0032] S7, globally search for the maximum value in the four-dimensional imaging data volume, where the maximum value represents the location of the microseismic event, and the time corresponding to the maximum value represents the occurrence time of the microseismic event. This achieves precise localization of microseismic events, providing crucial spatiotemporal distribution information for evaluating hydraulic fracturing effects, depicting fracture network morphology, and predicting production capacity. It solves the problems existing in traditional localization methods and improves the accuracy and efficiency of microseismic monitoring.

[0033] In this embodiment of the invention, the following steps are also included: On the CPU side, the three-dimensional medium model is decomposed into regions using a message passing interface, and the computational tasks are distributed to multiple GPU nodes. The microseismic data and the model parameters of the three-dimensional arbitrary curve physical mesh are stored and processed using a standard 32-bit single-precision floating-point format. On the GPU side, CUDA is used to perform the first-order velocity-stress equation solution and imaging condition calculation; In each thread, a local register of type FP32 is declared as an accumulator. When performing finite difference operations, the multiplication part is completed in BF16, and the multiplication-accumulation process is performed in this FP32 accumulator. After performing backpropagation iterations at each time step on the time axis, the reverse-time reconstructed wave field obtained at the next moment needs to be converted from the format of the BF16 or FP32 accumulator back to the FP32 format before being stored when writing the reverse-time reconstructed wave field back to the GPU global memory.

[0034] This invention, by introducing a BF16 mixed-precision computation strategy and CPU+GPU heterogeneous parallelism, can increase the computation speed of seismic source reverse-time imaging by several times while reducing memory and GPU memory usage by nearly half. This makes it possible to process larger-scale data and perform more detailed mesh simulations, significantly shortening the processing cycle and reducing hardware costs. This invention employs the BF16 format, supplemented by a high-precision accumulator and an FP32 master-replica storage strategy, successfully applying mixed-precision technology to solving complex wave equations. It enjoys the performance advantages of low-precision computation while ensuring numerical stability and physical reliability of the final results over long-term iterations. This invention integrates the two major challenges of efficient computation and low storage usage into a unified framework, forming a complete technical process that can be directly applied to production practice, greatly promoting the transition of seismic source reverse-time imaging technology from theoretical research to industrial application.

[0035] See Figure 2 Based on the above method, the present invention also discloses a microseismic reverse-time imaging system based on BF16 hybrid precision, comprising: A three-dimensional arbitrary curve physical mesh acquisition module is used to generate a three-dimensional arbitrary curve physical mesh based on surface elevation data and store the three-dimensional arbitrary curve physical mesh in FP32 format. The reverse-time reconstructed wavefield acquisition module is used to read the microseismic data recorded by the station in FP32 format, reverse the microseismic data on the time axis, and obtain the reverse-time reconstructed wavefield. The elastic wave propagation result acquisition module is used to reconstruct the wave field in reverse time as the boundary source term. On a three-dimensional arbitrary curve physical grid, the time iterative core calculation is performed on the GPU using CUDA in the BF16 format and the curve grid finite difference method to perform time iterative core calculations on all time steps. The first-order velocity-stress equations are solved iteratively to obtain the wave field master copy. The acceleration field decomposition module is used to decouple and separate the longitudinal and transverse waves in the main copy of the wave field in BF16 format to obtain the longitudinal wave acceleration field and the transverse wave acceleration field. The imaging value acquisition module at the grid point is used to acquire the BF16 format imaging value at the grid point based on the longitudinal wave acceleration field and the transverse wave acceleration field and the imaging condition formula. The four-dimensional imaging data volume acquisition module is used to aggregate the imaging values ​​at the grid points at all time steps, and transmit the aggregation result back to the CPU in FP32 format to construct the four-dimensional imaging data volume. The microseismic event confirmation module is used to globally search for the maximum value in the four-dimensional imaging data volume. The maximum value is the location of the microseismic event, and the time corresponding to the maximum value is the time when the microseismic event occurred.

[0036] The various modules of the system in this invention work together to achieve efficient computation and low storage usage for microseismic reverse time imaging.

[0037] Example 2: See Figure 3 The core objective of this invention is to provide a microseismic reverse-time imaging method based on BF16 hybrid precision, aiming to solve the problem of prominent computational resource bottlenecks in the background technology.

[0038] This invention employs a hybrid precision computing strategy based on a novel half-precision floating-point format, aiming to significantly reduce memory and video memory consumption during the reverse-time imaging process of seismic sources and significantly improve computing speed.

[0039] This invention designs a stable and reliable CPU+GPU heterogeneous parallel computing process through hardware and software co-optimization. While ensuring computing accuracy, it maximizes the utilization efficiency of computing resources, making large-scale, high-precision microseismic reverse-time imaging feasible for practical applications.

[0040] This invention first employs curved mesh technology to accurately simulate seismic wave propagation on undulating surfaces from a physical perspective; secondly, at the high-performance computing level, it introduces and optimizes the BF16 mixed-precision format to address bottlenecks in computational efficiency and resource consumption; finally, a sophisticated numerical implementation strategy ensures the stability and accuracy of low-precision calculations. The specific technical steps of this invention are as follows: Step 1: Generation of computational grid and model building including terrain; Input three-component microseismic records, surface elevation data (DEM), and subsurface P-wave velocity, S-wave velocity, and density models. Based on the surface elevation data, generate a three-dimensional arbitrary curve physical grid that can accurately fit the surface undulation, and establish the mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system.

[0041] Step 2, Wavefield Reverse Time Reconstruction Based on Curved Grid; Microseismic data recorded by the stations Reversing the time axis yields a reverse-time reconstructed wavefield. The process can be represented as:

[0042] in, From the position of the station to any point in space Green's function; This is the total duration of the data; Represents convolution; To reconstruct the wave field in reverse time; Spatial location coordinates; It is a time variable; Location of the station; Component data recorded by the station.

[0043] Using the reconstructed wave field in reverse time as the boundary source term, the elastic wave propagation described by the following first-order velocity-stress equations is solved in reverse time on the three-dimensional arbitrary curve physical mesh generated in step 1 using the high-order finite difference method:

[0044] in, Density; For speed Quantity; For speed Quantity; For speed Quantity; These are the components of the stress tensor; and For Lamé parameters; for Coordinates of direction; for Coordinates of direction; for The coordinates of the direction. By using the mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system, the above equations can be transformed into the curve grid coordinate system.

[0045] Step 3, Application of imaging conditions based on propagation direction constraints; In the wavefield reverse-time reconstruction calculation at each time step, for each grid point, based on the decoupled P-wave acceleration field... and transverse wave acceleration field Calculate the corresponding energy propagation direction vector. Energy propagation direction vector of the longitudinal wave acceleration field. The component form is:

[0046] Energy propagation direction vector of transverse wave acceleration field The component form is:

[0047] in, For the permutation tensor; This is the energy propagation direction vector of the longitudinal wave acceleration field; This is the energy propagation direction vector of the transverse wave acceleration field; Let be the divergence of the velocity field; For longitudinal wave acceleration wave field; and For the permutation tensor; For transverse wave acceleration wave field; Spatial derivative of the velocity field.

[0048] Subsequently, the energy propagation direction vector of the longitudinal wave acceleration field is applied. and the energy propagation direction vector of the transverse wave acceleration field The conditions for dot product imaging are mathematically expressed as follows:

[0049] in, For grid points At any moment The imaging values ​​are obtained by utilizing the physical property that the propagation directions of P-waves and S-waves are approximately parallel at the actual seismic source, which can effectively suppress imaging artifacts and improve imaging resolution.

[0050] Step 4: Hybrid precision heterogeneous parallel computing based on BF16; To efficiently execute steps 2 and 3, this invention employs a CPU+GPU heterogeneous parallel computing framework and designs an ingenious mixed-precision memory optimization strategy. Step 4.1, Task Division and Initialization: On the CPU side, the 3D medium model is decomposed into regions using the Message Passing Interface (MPI), and the computational tasks are distributed to multiple GPU nodes. During this stage, key variables such as microseismic data and model parameters of the 3D arbitrary curve physical mesh are stored and processed using standard 32-bit single-precision floating-point (FP32) format to ensure initial accuracy.

[0051] Step 4.2, GPU Core Computation and BF16 Optimization: On the GPU side, CUDA is used to execute the core loop for solving the first-order velocity-stress equations and calculating imaging conditions. To conserve valuable GPU memory and improve computational and memory access efficiency, wavefield variables are primarily calculated in 16-bit BF16 (Bfloat16) format on the GPU. The BF16 format has the same dynamic numerical range as FP32, effectively avoiding the numerical overflow problem that easily occurs in wavefield multiplication operations using the conventional half-precision (FP16) format.

[0052] Step 4.3, High-Precision Accumulator Strategy: To mitigate the accumulated errors that may arise from the relatively low mantissa precision of BF16, a high-precision accumulator strategy is implemented within the CUDA kernel functions. Specifically, an FP32-type local register is declared as the accumulator in each thread. When performing finite difference (weighted summation) operations, the multiplication part can be completed in BF16, but the multiplication and accumulation process is performed in this FP32 accumulator. This ensures numerical stability, trading significant precision for minimal performance overhead.

[0053] Step 4.4, Mixed-precision state storage: After each time step iteration, the wave field calculated for the next time step is reconstructed in reverse time. When writing back to the GPU global memory, it is converted from BF16 or FP32 accumulator back to FP32 format for storage. That is, the master copy of the wave field is always stored in global memory in high-precision FP32 format, thereby blocking the propagation and accumulation of errors between time steps.

[0054] Step 5: Image result aggregation and seismic source localization; After all time steps are calculated, a four-dimensional imaging data volume covering the entire space and time is formed by applying imaging conditions in step 3. By globally searching for the maximum value in this data volume, the corresponding spatial coordinates are the location of the microseismic event, and the corresponding time is the moment of occurrence of the microseismic event.

[0055] The present invention significantly improves computational efficiency and resource utilization by introducing a BF16 mixed-precision computing strategy and CPU+GPU heterogeneous parallelism. Figure 3 This method can increase the computation speed of seismic source reverse-time imaging by several times, while reducing memory and video memory usage by nearly half. Figure 8 This makes it possible to process larger-scale data and perform more detailed grid simulations, significantly shortening the processing cycle and reducing hardware costs. Figure 8 In the table, a represents the comparison of computation time, where the BF16 mixed-precision method takes 54.9% of the time of the FP32 precision method, and b represents the comparison of peak video memory usage, where the BF16 mixed-precision method uses 67.7% of the video memory of the FP32 precision method.

[0056] This invention exhibits high numerical stability and reliability: it selects the BF16 format and employs a high-precision accumulator and an FP32 primary replica storage strategy. Figure 4 This successfully applied mixed-precision techniques to solving complex wave equations, enjoying the performance advantages of low-precision computation. Figure 8 This ensures both numerical stability over long-term iterations and the physical reliability of the final result. Figure 5 , Figure 6 , Figure 7 ). Figure 6The results of microseismic reverse-time imaging localization tests on a 3D complex medium model are shown. a) is the source localization result based on BF16 hybrid accuracy reverse-time imaging, and b) is the source localization result based on conventional FP32 accuracy reverse-time imaging. The localization results from the two methods are basically consistent, and the locations coincide with the actual microseismic locations (red circles). Figure 7 The image shows the error histograms of the BF16 mixed-precision imaging results and the conventional FP32 precision imaging results for a three-dimensional complex medium model. The overall error follows a normal distribution, and most of the errors are within ±10%, which is relatively small.

[0057] This invention is highly practical and versatile: It integrates the two major challenges of efficient computing and low storage usage into a unified framework, forming a complete technical process that can be directly applied to production practice, which greatly promotes the transition of seismic source reverse time imaging technology from theoretical research to industrial application.

[0058] Example 3: See Figure 3 This embodiment provides a microseismic reverse-time imaging method based on BF16 hybrid precision, including the following steps: Step 1: Generation of computational grid and model building with terrain; First, publicly available terrain data for the work area is acquired and used as the upper boundary of the model. Combined with a one-dimensional layered velocity model, a three-dimensional medium model is generated, measuring 7km × 8km horizontally and covering the target depth vertically. Then, based on this three-dimensional medium model, a three-dimensional arbitrary curve mesh with a spatial grid spacing of 15m × 15m × 15m is generated, with the upper surface of this mesh precisely conforming to the actual terrain.

[0059] Step 2: Wavefield reverse propagation based on curved grid; 251 triadic microseismic data points from a single microseismic event were selected, with a total data duration of 5 seconds and a time sampling interval of 1 millisecond. The data from these 5000 time sampling points were time-reversed. After data reading and preprocessing on the CPU, the reversed data and curve mesh data were transferred to the memory of four GPUs. The reversed data was used as the source term and injected at the station location. A CUDA-based curve mesh finite difference solver was then used to perform inverse-time reconstruction of the wavefield.

[0060] Step 3: Application of Imaging Conditions Based on Propagation Direction Constraints In each iterative calculation on the GPU, in addition to updating the velocity and stress wave fields, the program simultaneously calculates the decoupled P-wave and S-wave acceleration fields, and further calculates their energy propagation direction vectors. and Then, the imaging value is calculated at each grid point. and the current time Imaging values With this grid point The stored historical maximum imaging values ​​are compared and updated.

[0061] Step 4: Implementation of mixed-precision heterogeneous parallel computing based on BF16; This step is the core of efficient calculation.

[0062] 4.1 Task Partitioning and Initialization: In the CPU main process, MPI is used to decompose the 3D computational mesh into four sub-regions along the Y-axis, which are then allocated to four GPUs. The initial velocity and density models and seismic records are loaded in FP32 format on the CPU and then distributed and transmitted to each GPU.

[0063] 4.2 GPU Core Computation: On each GPU, a master copy of the wavefield arrays, including velocity and stress, is stored in global memory in FP32 format. When a CUDA kernel function is launched to perform a time step calculation, the FP32 wavefield values ​​of the required neighboring points are read from global memory.

[0064] 4.3 Application of the high-precision accumulator strategy: The read FP32 values ​​are converted to BF16 format in the registers inside the kernel function. The multiplication part of the finite difference operation is performed using BF16 instructions, but the accumulation of all product terms is done in a pre-declared FP32 type local accumulator.

[0065] 4.4 Mixed-Precision State Storage: After the finite difference operation is completed, the result in the FP32 accumulator (i.e., the wavefield update for the next time step) is written back to the FP32 wavefield master copy in global memory. In this way, the computationally intensive parts utilize the high throughput of BF16, while the preservation and propagation of the state maintain the precision of FP32, effectively preventing error accumulation.

[0066] Step 5: Imaging result aggregation and seismic source localization After a complete reverse-time reconstruction over 5000 time steps, each GPU stored the maximum 3D imaging value data volume for its assigned region. MPI was used to reduce and stitch these four sub-data volumes back into CPU main memory, forming the complete final 3D imaging result. By searching for the global maximum point within this 3D data volume, the final location of the microseismic event was determined to be (X: 4.025km, Y: 3.877km, Z: -1.169km), with an occurrence time of 1.66 seconds.

[0067] Compared with conventional methods of pure FP32 computation, the embodiments of the present invention have achieved significant results, as detailed below: Accuracy Verification: Based on the positioning results of this invention, compared with the positioning results of conventional methods using pure FP32 calculations, the location of the seismic source is basically consistent with the time of earthquake occurrence. The statistical results of data errors for all imaging points in the imaging results are less than 10%, fully demonstrating that the method of this invention has the same imaging accuracy as the pure FP32 method.

[0068] Efficiency Improvement: Under the same hardware conditions, the computation time of this invention is approximately 2.5 hours, while the conventional FP32 method takes approximately 4.8 hours, representing a near doubling of computational efficiency. Simultaneously, peak memory usage is reduced from approximately 38GB to approximately 23.5GB, saving 39% of memory resources.

[0069] This invention can effectively solve the problem of accurate microseismic location in complex terrain, and can significantly improve computational efficiency and reduce resource consumption, thus having high practical application value.

[0070] Example 4: See Figure 1 and Figure 3 This embodiment provides a microseismic reverse-time imaging method based on BF16 mixed precision. This embodiment is an application example of microseismic monitoring data from hydraulic fracturing in shale gas extraction in a certain region. Figure 9 It is the velocity model of P-waves and S-waves in this region and the topographic relief. Figure 10 This is the input microseismic monitoring data for this embodiment. The signal-to-noise ratio of this microseismic data is low, but significant microseismic events can still be observed. Figure 11 The figures show the results of locating the microseismic data in this embodiment using the source reverse time imaging method based on BF16 mixed accuracy and the source reverse time imaging method based on pure FP32 accuracy proposed in this invention. As can be seen from the figures, the location results of the two methods are basically consistent. Figure 12 This is an error histogram comparing the imaging results using BF16 mixed precision and conventional pure FP32 precision imaging results in this embodiment. It can be seen that the overall errors of both follow a normal distribution, with most errors within ±10%, which is well within acceptable limits. Figure 13This example compares the computation time and video memory usage of the two methods. (a) shows the computation time comparison: the BF16 mixed-precision method takes 52.1% of the time of the FP32 precision method. (b) shows the peak video memory usage comparison: the BF16 mixed-precision method uses 61.8% of the video memory used by the FP32 precision method. Under the same hardware conditions, the computation time of this invention is approximately 2.5 hours, while the conventional FP32 method takes approximately 4.8 hours, nearly doubling the computational efficiency. Simultaneously, the peak video memory usage is reduced from approximately 38GB to approximately 23.5GB, consuming only 61.8% of the computer video memory resources of the conventional method. This demonstrates that the method of this invention can significantly improve computational efficiency and reduce resource consumption while maintaining accuracy, possessing high practical value.

[0071] An electronic device includes: a processor; a memory for storing computer program instructions; and for implementing a microseismic reverse-time imaging method based on BF16 mixed precision when executing the computer program.

[0072] A storage medium storing computer program instructions, which are loaded and executed by a processor, wherein the processor performs a microseismic reverse-time imaging method based on BF16 mixed precision.

[0073] A computer program product comprising computer instructions that instruct a computer to execute a microseismic reverse-time imaging method based on BF16 mixed precision.

[0074] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0075] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1A device that provides the functions specified in one or more boxes.

[0076] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0077] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0078] The above content is only for illustrating the technical concept of the present invention and should not be construed as limiting the scope of protection of the present invention. Any modifications made to the technical solution based on the technical concept proposed in this invention shall fall within the scope of protection of this invention.

Claims

1. A microseismic reverse-time imaging method based on BF16 hybrid precision, characterized in that, Using a hybrid BF16 and FP32 storage format, and employing a heterogeneous parallel computing framework consisting of CPUs and GPUs, the following steps are included: A three-dimensional arbitrary curve physical mesh is generated based on surface elevation data, and the three-dimensional arbitrary curve physical mesh is stored in FP32 format. The microseismic data recorded by the station was read in FP32 format, and the microseismic data was reversed on the time axis to obtain the reverse time reconstructed wavefield; Using the reconstructed wavefield in reverse time as the boundary source term, on a three-dimensional arbitrary curved physical mesh, the finite difference method of the curved mesh in BF16 format is used to perform time iteration core calculations on all time steps on the GPU through CUDA, and the first-order velocity-stress equations are solved iteratively to obtain the wavefield master copy. The P-wave and S-wave acceleration fields are decoupled and separated from the primary and secondary copies of the wave field using the BF16 scheme to obtain the P-wave acceleration field and the S-wave acceleration field. Based on the longitudinal wave acceleration field and the transverse wave acceleration field, and using the imaging condition formula, the BF16 format imaging values ​​at the grid points are obtained. The imaging values ​​at the grid points of all time steps are aggregated, and the aggregation result is transmitted back to the CPU in FP32 format to construct a four-dimensional imaging data volume. The maximum value in the four-dimensional imaging data volume is searched globally. The maximum value is the location of the microseismic event, and the time corresponding to the maximum value is the time when the microseismic event occurred.

2. The microseismic reverse-time imaging method based on BF16 hybrid precision according to claim 1, characterized in that, The formula for obtaining the reconstructed wave field in reverse time is as follows: in, From the position of the station to any point in space Green's function; This is the total duration of the data; Represents convolution; To reconstruct the wave field in reverse time; Spatial location coordinates; It is a time variable; Location of the station; Component data recorded by the station.

3. The microseismic reverse-time imaging method based on BF16 hybrid precision according to claim 1, characterized in that, The process of performing time-iteration core calculations at all time steps, iteratively solving the first-order velocity-stress equations to obtain the wavefield master copy includes: A mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system is established based on the aforementioned three-dimensional arbitrary curve physical mesh; By using the mapping relationship between the physical space coordinate system and the computational space Cartesian coordinate system, the first-order velocity-stress equations are transformed into the curve grid coordinate system to obtain the transformed first-order velocity-stress equations. Based on the transformed first-order velocity-stress equations, at each time step, the wave field value is read from the FP32 global memory, the wave field value is converted to BF16 format in the FP32 local accumulator, multiplication is performed and the value is accumulated to obtain the accumulator result, and then the accumulator result is converted back to FP32 format and written back to the global memory as the main copy of the wave field.

4. The microseismic reverse-time imaging method based on BF16 hybrid precision according to claim 1, characterized in that, The first-order velocity-stress equations are as follows: in, Density; For speed Quantity; For speed Quantity; For speed Quantity; These are the components of the stress tensor; and For Lamé parameters; for Coordinates of direction; for Coordinates of direction; for Coordinates of direction.

5. The microseismic reverse-time imaging method based on BF16 hybrid precision according to claim 3, characterized in that, It also includes the following steps: On the CPU side, the three-dimensional medium model is decomposed into regions using a message passing interface, and the computational tasks are distributed to multiple GPU nodes. The microseismic data and the model parameters of the three-dimensional arbitrary curve physical mesh are stored and processed using a standard 32-bit single-precision floating-point format. On the GPU side, CUDA is used to perform the first-order velocity-stress equation solution and imaging condition calculation; In each thread, a local register of type FP32 is declared as an accumulator. When performing finite difference operations, the multiplication part is completed in BF16, and the multiplication-accumulation process is performed in this FP32 accumulator. After performing backpropagation iterations at each time step on the time axis, the reverse-time reconstructed wave field obtained at the next moment needs to be converted from the format of the BF16 or FP32 accumulator back to the FP32 format before being stored when writing the reverse-time reconstructed wave field back to the GPU global memory.

6. The microseismic reverse-time imaging method based on BF16 hybrid precision according to claim 1, characterized in that, The step of obtaining the imaging values ​​at the grid points based on the P-wave acceleration field and the S-wave acceleration field and the imaging condition formula includes: The energy propagation direction vector of the longitudinal wave acceleration field is obtained from the longitudinal wave acceleration field. Obtain the energy propagation direction vector of the shear wave acceleration field from the shear wave acceleration field; Based on the energy propagation direction vectors of the longitudinal wave acceleration field and the transverse wave acceleration field, and using the imaging condition formula, the imaging values ​​at the grid points are obtained. The imaging condition formula is as follows: in, For grid points At any moment The imaging value; This is the energy propagation direction vector of the longitudinal wave acceleration field; This is the energy propagation direction vector of the transverse wave acceleration field.

7. The microseismic reverse-time imaging method based on BF16 hybrid precision according to claim 6, characterized in that, The formula for obtaining the energy propagation direction vector of the longitudinal wave acceleration field is: The formula for obtaining the energy propagation direction vector of the transverse wave acceleration field is: in, For the permutation tensor; This is the energy propagation direction vector of the longitudinal wave acceleration field; This is the energy propagation direction vector of the transverse wave acceleration field; Let be the divergence of the velocity field; For longitudinal wave acceleration wave field; and For the permutation tensor; For transverse wave acceleration wave field; Spatial derivative of the velocity field.

8. A microseismic reverse-time imaging system based on BF16 hybrid precision, characterized in that, include: A three-dimensional arbitrary curve physical mesh acquisition module is used to generate a three-dimensional arbitrary curve physical mesh based on surface elevation data and store the three-dimensional arbitrary curve physical mesh in FP32 format. The reverse-time reconstructed wavefield acquisition module is used to read the microseismic data recorded by the station in FP32 format, reverse the microseismic data on the time axis, and obtain the reverse-time reconstructed wavefield. The elastic wave propagation result acquisition module is used to reconstruct the wave field in reverse time as the boundary source term. On a three-dimensional arbitrary curve physical grid, the time iterative core calculation is performed on the GPU using CUDA in the BF16 format and the curve grid finite difference method to perform time iterative core calculations on all time steps. The first-order velocity-stress equations are solved iteratively to obtain the wave field master copy. The acceleration field decomposition module is used to decouple and separate the longitudinal and transverse waves in the main copy of the wave field in BF16 format to obtain the longitudinal wave acceleration field and the transverse wave acceleration field. The imaging value acquisition module at the grid point is used to acquire the BF16 format imaging value at the grid point based on the longitudinal wave acceleration field and the transverse wave acceleration field and the imaging condition formula. The four-dimensional imaging data volume acquisition module is used to aggregate the imaging values ​​at the grid points at all time steps, and transmit the aggregation result back to the CPU in FP32 format to construct the four-dimensional imaging data volume. The microseismic event confirmation module is used to globally search for the maximum value in the four-dimensional imaging data volume. The maximum value is the location of the microseismic event, and the time corresponding to the maximum value is the time when the microseismic event occurred.

9. An electronic device, comprising: A processor; a memory, an electronic device for storing computer program instructions; characterized in that, when executing the computer program, it implements the microseismic reverse-time imaging method based on BF16 hybrid precision as described in any one of claims 1-7.

10. A storage medium storing computer program instructions, characterized in that, When the computer program instructions are loaded and run by the processor, the processor executes the microseismic reverse-time imaging method based on BF16 mixed precision as described in any one of claims 1-7.

Images