A GPU-based method, device, and medium for rapid inversion of 3D ore bodies.

Abstract

This invention discloses a GPU-based method, device, and medium for rapid 3D orebody inversion, relating to the field of 3D orebody inversion technology. The method includes importing multi-source exploration data, generating observation data in a unified coordinate system, establishing a 3D voxel grid, and obtaining a data acquisition configuration package; performing multiphysics forward modeling to generate normalized prediction vectors and normalized residual vectors, establishing the forward and transpose interfaces of the sensitivity operator, constructing update direction vectors, and generating a pre-projection parameter field; sequentially performing box-constrained projection, total variation projection, and wavelet sparse projection on the pre-projection parameter field to obtain a sparse parameter field and generate the final parameter field; generating a 3D grade voxel array, calculating the grade segmentation threshold, extracting the orebody connected body set, generating a orebody triangular grid and an orebody existence probability voxel array, and obtaining a digital twin base data package. This method shortens computation time, reduces storage pressure, and improves the accuracy and practicality of the inversion.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional ore body inversion technology, and in particular to a GPU-based method, equipment, and medium for rapid three-dimensional ore body inversion. Background Technology

[0002] The core objective of 3D orebody inversion is to infer the spatial morphology and grade distribution of underground orebodies by integrating multi-source observation data such as magnetic, electrical, induced polarization, seismic, and well logging data. The routine process includes importing and registering multi-source data and coordinates, characterizing and analyzing observation errors, constructing a physical property parameter field on a 3D voxel grid, and using multiphysics forward modeling and sensitivity operator-driven iterative updates to ensure that the predicted response and observation vectors meet statistical consistency, thereby inverting the spatial structure and grade information of the orebody.

[0003] However, the conventional process still has room for improvement in two aspects. First, multiphysics forward modeling and the calculation of associated or transpose effects often account for the main computational overhead. As the size of the 3D mesh increases, the iteration time of forward modeling and gradient calculation increases sharply, and the pressure on memory and storage increases significantly. If the borehole property constraints and spatial priors lack a unified and coordinated projection processing mechanism, it is easy to cause an imbalance in the application of various constraints, making it difficult to control the iterative convergence process and affecting the stability and reliability of the inversion. Summary of the Invention

[0004] In view of the aforementioned existing problems, the present invention is proposed.

[0005] Therefore, this invention provides a GPU-based method for fast inversion of 3D ore bodies to solve the problems of time-consuming physical forward iterative calculations, high storage consumption, and difficulty in unifying the projection fusion of constraint priors in existing technologies.

[0006] To solve the above-mentioned technical problems, the present invention provides the following technical solution:

[0007] In a first aspect, the present invention provides a GPU-based method for rapid inversion of 3D ore bodies, comprising: importing multi-source exploration data, generating observation data in a unified coordinate system, establishing a 3D voxel grid, generating a borehole index array and a box-constrained interval table of borehole properties, generating an initial parameter field, and obtaining an acquisition configuration package; based on the 3D voxel grid and the acquisition configuration package, performing multiphysics forward modeling, generating a normalized prediction vector and a normalized residual vector, establishing a forward action interface and a transpose action interface of the sensitivity operator, constructing an update direction vector, and generating a pre-projection parameter field; generating a trigger set based on the acquisition configuration package, sequentially performing box-constrained projection, total variation projection, and wavelet sparse projection on the pre-projection parameter field to obtain a sparse parameter field, performing box-constrained projection based on the sparse parameter field to generate a final parameter field; generating a 3D grade voxel array based on the final parameter field, calculating the grade segmentation threshold, extracting the ore body connected body set, generating an ore body triangular grid and an ore body existence probability voxel array, and obtaining a digital twin base data package.

[0008] As a preferred embodiment of the GPU-based rapid inversion method for 3D ore bodies described in this invention, the specific steps for obtaining the acquisition configuration package are as follows: importing multi-source exploration data, calculating coordinate transformation parameters, obtaining observation data under a unified coordinate system, and calculating the standard deviation vector for each type of data; based on the observation data under the unified coordinate system, calculating the horizontal and vertical grid step sizes to generate a 3D voxel grid; in the 3D voxel grid, searching for the nearest voxel for each borehole trajectory point to generate a borehole index array, and combining it with the standard deviation vector to generate a box-constrained interval table of borehole properties; based on the borehole index array and the box-constrained interval table of borehole properties, interpolating on the 3D voxel grid using Gaussian radial basis functions to generate an initial parameter field and obtain the acquisition configuration package.

[0009] As a preferred embodiment of the GPU-based fast inversion method for 3D ore bodies described in this invention, the specific steps for generating normalized prediction vectors and normalized residual vectors are as follows: based on a 3D voxel grid and a data acquisition configuration package, perform multiphysics forward modeling to generate normalized prediction vectors for various types of data; subtract the normalized observation vectors of various types of data from the normalized prediction vectors element by element to obtain the normalized residual vectors.

[0010] As a preferred embodiment of the GPU-based fast inversion method for 3D ore bodies described in this invention, the generation of the pre-projection parameter field includes: establishing a positive action interface and a transpose action interface for the sensitivity operator pair vector; calculating the checkpoint interval; calling the transpose action interface based on the normalized residual vector; constructing the update direction vector using a rotating master-driven approach; calling the positive action interface; calculating the parameter update step size; and performing an update along the update direction vector based on the current iterative parameter field to generate the pre-projection parameter field.

[0011] As a preferred embodiment of the GPU-based rapid inversion method for three-dimensional ore bodies described in this invention, the generation of the trigger set includes: constructing a set of physical observation types based on the acquisition configuration package, calculating the data consistency target norm, and generating the trigger set.

[0012] As a preferred embodiment of the GPU-based fast inversion method for 3D ore bodies described in this invention, the specific steps for generating the final parameter field are as follows: Based on the pre-projection parameter field, combined with the borehole index array and the box-constrained interval table of borehole properties, the voxel values ​​of the pre-projection parameter field at the borehole index positions are truncated in parallel to obtain the box-constrained parameter field; the residual ratio is calculated to obtain the total variation projection target scalar, and the dual variable iteration is performed on the box-constrained parameter field to generate the total variation parameter field; a unified noise scale is calculated based on the trigger set, a soft thresholding threshold is calculated, a 3D orthogonal wavelet transform is performed on the total variation parameter field and soft thresholding is performed, and then an inverse wavelet transform is performed to generate a sparse parameter field; based on the sparse parameter field, box-constrained projection is performed through the borehole index array and the box-constrained interval table of borehole properties to generate the final parameter field.

[0013] As a preferred embodiment of the GPU-based rapid inversion method for three-dimensional ore bodies described in this invention, the specific steps for extracting the ore body connected component set are as follows: based on the final parameter field, the borehole index array is scanned in parallel to construct a physical property grade training sample table; a grade prediction operator is established using kernel regression to generate a three-dimensional grade voxel array; the grade segmentation threshold is calculated to generate a binary ore body indicator voxel array; three-dimensional connected component labeling is performed to extract the ore body connected component set.

[0014] As a preferred embodiment of the GPU-based rapid inversion method for 3D ore bodies described in this invention, the specific steps for obtaining the digital twin base data package are as follows: for each ore body in the ore body connectivity set, extract the boundary voxels of the ore body connectivity to generate an ore body triangular mesh; construct a perturbation observation vector set based on the standard deviation vector, perform a single round of outer layer iteration to obtain an ore body existence probability voxel array, and combine the 3D grade voxel array and the ore body triangular mesh to obtain the digital twin base data package.

[0015] In a second aspect, the present invention provides a computer device including a memory and a processor, wherein the memory stores a computer program, wherein when the computer program is executed by the processor, it implements any step of the GPU-based fast inversion method for three-dimensional ore bodies as described in the first aspect of the present invention.

[0016] Thirdly, the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein: when the computer program is executed by a processor, it implements any step of the GPU-based fast inversion method for three-dimensional ore bodies as described in the first aspect of the present invention.

[0017] The beneficial effects of this invention are as follows: by using a GPU parallel computing framework to process multiple data blocks simultaneously, the speed of large-scale data processing in the three-dimensional ore body inversion process is improved, the calculation time is shortened, the storage pressure is reduced, and the processing efficiency is improved; by using an optimization algorithm based on projection constraints to accurately fuse prior information and observation data, and by adaptively determining the soft thresholding threshold based on the trigger set and unified noise scale, the convergence controllability and inversion stability are improved. Attached Figure Description

[0018] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a flowchart of a GPU-based method for rapid inversion of 3D ore bodies.

[0020] Figure 2 This is a graph showing the average time taken for a single outer layer iteration versus the total number of voxels.

[0021] Figure 3 The graph shows the normalized residual vector L2 norm and outer iteration number for each observation type.

[0022] Figure 4 A flowchart for obtaining the data acquisition configuration package.

[0023] Figure 5 A flowchart for generating the parameter field before projection.

[0024] Figure 6 A flowchart for obtaining the final parameter field. Detailed Implementation

[0025] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0026] Many specific details are set forth in the following description in order to provide a full understanding of the invention. However, the invention may also be practiced in other ways different from those described herein, and those skilled in the art can make similar extensions without departing from the spirit of the invention. Therefore, the invention is not limited to the specific embodiments disclosed below.

[0027] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The phrase "in one embodiment" appearing in different places in this specification does not necessarily refer to the same embodiment, nor is it a single or selective embodiment that is mutually exclusive with other embodiments.

[0028] Reference Figures 1-6 This is one embodiment of the present invention, which provides a GPU-based method for fast inversion of three-dimensional ore bodies, comprising the following steps:

[0029] S1. Import multi-source exploration data, generate observation data under a unified coordinate system, establish a three-dimensional voxel grid, generate a borehole index array and a box-constrained interval table of borehole properties, generate an initial parameter field, and obtain the acquisition configuration package.

[0030] Import multi-source exploration data, calculate coordinate transformation parameters, obtain observation data under a unified coordinate system, and calculate the standard deviation vector of each type of data; based on the observation data under the unified coordinate system, calculate the horizontal and vertical grid step size, and generate a three-dimensional voxel grid.

[0031] Furthermore, magnetic data, electrical data, induced polarization data, seismic data, well logging data, and borehole data are imported using independent table structures. Each type of data table includes, but is not limited to, observation value fields, observation point coordinate fields, observation type identifier fields, and acquisition batch identifier fields. The import process adopts fixed parsing rules. Specifically, the semantics of the fields are first matched by column names, and then the fields are parsed and verified in sequence. Specifically, the coordinate field is parsed into three-dimensional coordinates, the observation value field is parsed into real numbers, and the acquisition batch identifier is parsed into a string and participates in the versioned metadata hash.

[0032] Missing fields discovered during the parsing process are marked with fixed missing tags. The missing field name, corresponding data record, and missing tag are written into the versioned metadata. Outliers in the observed value fields are marked but not removed. A dual-condition trigger mechanism is used to identify outliers. Specifically, when either the absolute amplitude out-of-bounds tag or the quantile range out-of-bounds tag is triggered, the observed value is marked as an outlier. The outlier tag information is associated with the index of the corresponding data record and saved synchronously with the versioned metadata.

[0033] Using the mining area survey benchmark as the coordinate frame, coordinate transformation parameters are obtained by rigid body registration with the benchmark points. Specifically, three or more pairs of benchmark points are read from the survey file to form a source point set and a target point set. The centroids of the source point set and the target point set are calculated and decentralized. Covariance matrix decomposition is performed on the decentralized point set, and the rotation matrix is ​​obtained through singular value decomposition. The translation vector is obtained through the centroid relationship. The rotation matrix and the translation vector are assembled into a homogeneous transformation matrix and written into versioned metadata. Coordinate transformation is performed on each magnetic observation point, electrical electrode point, induced polarization electrode point, seismic shot point and receiver point, and borehole trajectory point to generate unified coordinates. The minimum envelope box with unified coordinates is used as the mining area envelope and written into the versioned metadata.

[0034] The standard deviation vector is calculated, specifically, by including instrument error and repeated measurement spread. Instrument error is obtained by reading the acquisition records or calibration files; repeated measurement spread is obtained by calculating the sample standard deviation of repeated observations at the same measurement point. For each observation point, the instrument error and repeated measurement spread are combined by taking the square root of the sum of their squares to obtain the observation point standard deviation. The standard deviations of all observation points are written into a diagonal matrix and added to the acquisition configuration package. The method for generating the normalized observation vector is also written into the acquisition configuration package. For each type of observation vector, the normalized observation vector is obtained by dividing each element by the observation point standard deviation.

[0035] A three-dimensional voxel mesh is established. Specifically, in the horizontal plane, all observation points are projected onto the plane, and the horizontal distance from each observation point to its nearest neighbor is calculated. The median of the nearest neighbor horizontal distances for all observation points is used as the horizontal sampling interval. In the vertical direction, the absolute values ​​of adjacent differences in elevation or depth of observation points are calculated, and the median of the absolute values ​​of adjacent differences is used as the vertical sampling interval. Half of the horizontal sampling interval is selected as the horizontal step size, and half of the vertical sampling interval is selected as the vertical step size.

[0036] Based on the mining area envelope and grid step size, the grid dimension is calculated and rounded up. A voxel center coordinate array is generated in a fixed index order, and a mapping from 3D index to linear index is generated. The grid origin, grid step size, grid dimension, and index order are written into the versioned metadata.

[0037] Furthermore, in the three-dimensional voxel mesh, the nearest voxel is searched for each borehole trajectory point to generate a borehole index array. Combined with the standard deviation vector, a box-constrained interval table of borehole properties is generated. Based on the borehole index array and the box-constrained interval table of borehole properties, Gaussian radial basis functions are used to interpolate on the three-dimensional voxel mesh to generate an initial parameter field and obtain the acquisition configuration package.

[0038] A borehole index array is generated. Specifically, for each borehole trajectory point, using the grid origin and grid step size, the difference between the coordinates of the borehole trajectory point and the coordinates of the grid origin is divided by the grid step size and rounded to the nearest integer. The coarse localization voxel 3D index of the borehole trajectory point is calculated. A 3×3×3 cubic neighborhood is constructed with the coarse localization voxel as the center. The index offset of the cubic neighborhood in each direction is -1, 0, or +1, and it includes a maximum of 27 candidate voxels. The Euclidean distance from the borehole trajectory point to the center coordinates of each candidate voxel is calculated, and the candidate voxel with the smallest Euclidean distance is used as the borehole index of the borehole trajectory point. If the borehole trajectory point is located near the grid boundary, causing some candidate voxel indices to exceed the boundary, the candidate voxels that exceed the boundary are ignored. The selection is compared and matched within the valid candidate voxel range, and all borehole trajectory points are matched to generate the borehole index array.

[0039] Furthermore, based on the borehole index array, the logging curve or borehole property measurement value corresponding to each borehole index position is read. The property includes magnetic susceptibility, electrical conductivity, polarization parameters, and velocity. For each type of property, the standard deviation of the property is calculated according to the instrument accuracy or repeated measurement scatter information in the measurement file. A box-constrained interval table is generated with the property measurement value as the center and three times the standard deviation as the half-width. The lower bound is the difference between the measurement value and three times the standard deviation, and the upper bound is the sum of the measurement value and three times the standard deviation. The borehole index array and the box-constrained interval table of the borehole property are stored in the same index order and written into the acquisition configuration package. At the same time, the summary information is written into the versioned metadata.

[0040] It should be noted that the box constraint interval formed by three times the standard deviation and half the width has sufficient coverage for conventional noise, and three times the standard deviation can reduce the truncation outside the interval caused by occasional fluctuations, local outliers or transient sampling disturbances.

[0041] Furthermore, interpolation is performed using Gaussian radial basis functions. For each voxel center coordinate, the K nearest borehole index points are selected from all borehole index points to form a neighbor sample set. The value of K is determined by taking the square root of the total number of borehole index points and rounding it up.

[0042] The Gaussian radial basis function used is of exponential decay form, and the shape parameter is determined by the spatial distribution of borehole index points. Specifically, the nearest neighbor distance between all pairs of borehole index points is calculated, and twice the median of the nearest neighbor distance is used as the value of the shape parameter. The distance from the center of each voxel to the K neighboring borehole index points is calculated, and the original weight of each borehole index point is calculated using the Gaussian radial basis function. The original weights are normalized so that the sum of all weights is 1. The normalized weight of each neighboring borehole index point is multiplied by the corresponding borehole physical property measurement, and all products are summed to obtain the initial parameter value of the current voxel position. All voxels are traversed to generate the initial parameter field.

[0043] Specifically, the Gaussian radial basis function is expressed as:

[0044] ;

[0045] in, The Euclidean distance from the voxel center to the borehole sample point is denoted as . The original weights at that time This represents the Euclidean distance from the voxel center to the borehole sample point. Indicates shape parameters.

[0046] The acquisition configuration package includes, but is not limited to, magnetic, electrical, induced polarization (IP) and seismic, logging and borehole sections, as well as an error characterization section. Specifically, the magnetic section includes, but is not limited to, the dominant field direction, dominant field intensity, and observation type identifier; the electrical section includes, but is not limited to, electrode combinations and electrode coordinate indices, injection current identifiers, boundary condition types, and boundary parameters; the IPV section includes, but is not limited to, time window sets, time window start and end times, time window integration step size, IPV response kernel type, and kernel parameters; the seismic section includes, but is not limited to, sampling interval, record length, source wavelet description, and absorption boundary parameter set; the logging and borehole section includes, but is not limited to, borehole trajectory point sets, borehole index arrays, and box-constrained interval tables of borehole properties; and the error characterization section includes, but is not limited to, various standard deviation vectors and normalization calibers. The acquisition configuration package is serialized using fixed field names and field order, and the serialized text participates in the hash calculation of versioned metadata. The acquisition configuration package further includes an index table for storing observation type identifiers in the data space vector storage location and an index table for storing observation type identifiers in the standard deviation vector storage location. The index tables use the observation type identifier string as the key, and the key arrangement is fixed according to the serialized field order of the acquisition configuration package.

[0047] It should be noted that the voxel values ​​of magnetic susceptibility, conductivity, polarization parameters, and velocity are stored continuously using a structured array separation method; while the voxel center coordinate array, index mapping array, borehole index array, and box constraint interval table are stored as continuous arrays.

[0048] The magnetic susceptibility parameter field, electrical conductivity parameter field, polarization parameter field, and velocity parameter field on the three-dimensional voxel mesh are expanded and concatenated according to the voxel linear index order to obtain the parameter space vector. The concatenation order is the voxel value vector of magnetic susceptibility, the voxel value vector of electrical conductivity, the voxel value vector of polarization parameter, and the voxel value vector of velocity. After concatenation, a one-dimensional array is formed and used as the parameter space vector.

[0049] After sorting a certain type of physical observation according to the sorting keywords fixed in the acquisition configuration package, a one-dimensional array is generated to obtain a data space vector. The sorting keywords adopt a four-level ascending order, including the acquisition batch identifier, survey line or shot number, electrode combination number or point number, and time window number or time sampling number. Corresponding data space vectors are formed for magnetic, electrical, induced polarization, seismic, well logging, and borehole physical properties. For each type of data space vector, the starting offset, length, and sorting keyword caliber of the one-dimensional array are written into the index table of the acquisition configuration package, and the index table summary is written into the versioned metadata.

[0050] It should be noted that the parameter field is a set of physical property parameter fields discretely indexed by voxels on a three-dimensional voxel mesh. The parameter field includes the magnetic susceptibility parameter field, the electrical conductivity parameter field, the polarization parameter field, and the velocity parameter field. The initial parameter field is written to the versioned metadata iteration log field and stored as the parameter field entry for the current iteration. When the outer iteration first enters the multiphysics forward modeling and parameter update stage, it reads the parameter field entry for the current iteration from the versioned metadata iteration log field as the current iteration parameter field, and initializes the iteration number to zero.

[0051] Furthermore, the initial parameter field is copied element by element on the graphics processor according to the voxel linear index to generate a reference parameter field, and the reference parameter field is written to the reference field of the versioned metadata iteration log field.

[0052] Furthermore, based on the reference parameter field, the reference total variation scalar is calculated and written into the reference field of the versioned metadata iteration log field. Specifically, the three-dimensional differential gradient is calculated for the four voxel value arrays of the reference parameter field, the gradient magnitude is calculated and accumulated over all voxels, and the accumulated values ​​of the four physical property components are summed to obtain the reference total variation scalar. The three-dimensional differential gradient calculation uses a scale-consistent differential division with a three-dimensional voxel mesh step size.

[0053] Furthermore, to support the calculation of grade spatial distribution, the grade field is parsed in the borehole data table and written into the acquisition configuration package. Specifically, the test grade field is read from the borehole data table and parsed into a real number; the top depth field and bottom depth field of the test sample section are read from the borehole data table and parsed into real numbers; for each test sample section, the depth range is determined by the top depth field and bottom depth field of the test sample section; for each borehole trajectory point, the trajectory point depth is determined by the borehole trajectory point depth field; by determining whether the trajectory point depth falls within the depth range, the test grade field is associated with the borehole trajectory point; the associated grade value is written into the borehole grade array according to the index order of the borehole index array, and the borehole grade array is written into the logging and borehole section of the acquisition configuration package.

[0054] Furthermore, to support sparse projection, the three-dimensional orthogonal wavelet basis type field is written into the acquisition configuration package. The three-dimensional orthogonal wavelet basis type field is obtained by writing the string "Haar" before serializing the acquisition configuration package, and participates in the versioned metadata hash along with the serialized text of the acquisition configuration package.

[0055] S2. Based on the three-dimensional voxel mesh and acquisition configuration package, perform multiphysics forward modeling, generate normalized prediction vectors and normalized residual vectors, establish the forward action interface and transpose action interface of the sensitivity operator, construct the update direction vector, and generate the pre-projection parameter field.

[0056] Based on the 3D voxel mesh and the acquisition configuration package, multiphysics forward modeling is performed to generate normalized prediction vectors for various types of data. The normalized observation vectors of various types of data are subtracted element-wise from the normalized prediction vectors to obtain the normalized residual vectors.

[0057] Furthermore, based on the three-dimensional voxel mesh and the acquisition configuration package, a multi-physics forward modeling process is executed sequentially. Specifically, for magnetic forward modeling, a voxel volume discretization method is used, and the prediction vector is obtained through parallel reduction calculation of observation points on a graphics processor. For electrical forward modeling, a sparse linear equation system is formed by finite volume discretization, and the background electric potential field is solved using a conjugate gradient iterative algorithm on a graphics processor. The prediction vector is obtained by sampling the observation locations. For induced polarization forward modeling, the time discrete points are cyclically processed within each time window to construct the equivalent conductivity, and the electrical forward modeling solution process is repeated. The prediction vector is generated by averaging within the time window. For seismic forward modeling, a finite difference method is used for time advancement to simulate and obtain synthetic seismic records. The prediction vector is formed by sampling the locations of receiver points. For well logging and borehole property forward modeling, the corresponding voxel array is sampled through the borehole index to obtain the prediction vector. Each type of prediction vector is divided element-wise according to the standard deviation vector to generate a normalized prediction vector.

[0058] For each type of data, the normalized observation vector is subtracted element-wise from the normalized prediction vector to obtain the normalized residual vector, which is directly stored in the graphics processor memory. The L2 norm and infinity norm of each type of normalized residual vector are calculated using a reduction kernel function and written to the versioned metadata iteration log field. The L2 norm and infinity norm calculated by the reduction kernel function are written to the versioned metadata iteration log field with the observation type identifier string as the key, and the order of the keys is consistent with the key order of the index table from the observation type identifier of the acquisition configuration package to the data space vector storage location.

[0059] Furthermore, a positive action interface and a transpose action interface for the sensitivity operator on the vector are established, the checkpoint interval is calculated, the transpose action interface is called based on the normalized residual vector, the update direction vector is constructed using a rotating master-driven approach, the positive action interface is called, the parameter update step size is calculated, and an update is performed once along the update direction vector based on the current iteration parameter field to generate the parameter field before projection.

[0060] It should be noted that when the outer iteration enters the multi-physics forward modeling and parameter update stage, the current iteration parameter field and iteration number are read from the versioned metadata iteration log field; the iteration number is initialized to zero on the first entry, and the iteration number written last time is used on subsequent entries.

[0061] The sensitivity operator has a forward action interface and a transpose action interface for vectors. The interface includes a forward action interface and a transpose action interface. The forward action interface takes a parameter space vector as input and outputs a data space vector as output, and is used to calculate the linearized prediction increment. The transpose action interface takes a data space vector as input and outputs a parameter space vector as output, and is used to calculate the gradient and normal equation multipliers.

[0062] The interface implementation method is determined according to the physical type. Specifically, for seismic physics, the forward action interface is implemented using the first-order scattering approximation, which injects the increment of the input velocity parameter into the wave equation discretization to form a scattering source and performs a linearized forward modeling, outputting to the receiver point for sampling. The transpose action interface is implemented through adjoint injection, which injects the data space vector as the receiver point time series source into the adjoint equation and accumulates it by inner product with the forward wave field time series to obtain the parameter space vector of the velocity component.

[0063] For electrical methods and induced polarization physics, the forward action interface is implemented by solving for the first-order background potential and the incremental potential. The incremental potential is constructed from the equivalent source term caused by the increment of the input conductivity parameter or polarization parameter, and the sparse linear equation solver of the electrical forward model is reused. The equivalent source term is constructed by expressing the first-order perturbation of the increment of the conductivity parameter or polarization parameter using a discrete operator. The increment of the conductivity parameter or polarization parameter is combined with the background potential gradient on a voxel to generate the right-hand side vector used to solve for the incremental potential. The transpose action interface is implemented by solving for the first-order adjoint potential and accumulating it by performing a discrete inner product with the background gradient field to form a parameter space vector.

[0064] For magnetic physics, the forward-acting interface uses the increment of the input magnetic susceptibility parameter as the voxel weight and accumulates it according to the kernel function to obtain the observation increment; the transpose-acting interface back-projects the data space vector onto the voxel according to the kernel function to obtain the parameter space vector of the magnetic susceptibility component.

[0065] It should be noted that all interfaces correspond one-to-one with the fields in the acquisition configuration package. The fields in the acquisition configuration package that correspond one-to-one with the interfaces include, but are not limited to, the seismic sampling interval, record length, source wavelet and absorption boundary parameter set, electrode combination for electrical and induced polarization methods, injection current identifier, time window set and time window integration step size, magnetic field direction, main field intensity, and observation type identifier.

[0066] To support seismic wavefield co-computation, a checkpoint strategy is implemented to determine the storage interval of wavefield slices. The checkpoint interval is calculated by the available video memory of the graphics processor and the size of a single wavefield slice. Specifically, the number of bytes per wavefield slice is calculated based on the number of grid points and wavefield variables in the three-dimensional voxel. The wavefield variables are taken as two types in acoustic discretization: wavefield and wavefield first-order time derivative. The number of single-precision bytes is fixed at 4, and the number of slice bytes is the product of the number of grid points, the number of variables, and the number of single-precision bytes.

[0067] The available video memory of the graphics processor is dynamically read at runtime to avoid the video memory being fully occupied. The available video memory budget for the checkpoint is 35% of the available video memory of the graphics processor. The number of storable slices is obtained by using the ratio of the checkpoint budget to the number of bytes of a single slice and rounding down.

[0068] It should be noted that when repeatedly running 3D inversion on the same graphics processor, the runtime memory usage includes checkpoints, wave field propagation cache, gradient accumulation buffer, convolution and interpolation temporary regions, instantaneous allocation of parallel kernel functions, and memory fragmentation. If the checkpoint budget is too high, the peak superposition of temporary regions in the later stages of iteration can easily trigger overflow; if the checkpoint budget is too low, the insufficient number of slices will lead to frequent recalculations and a significant increase in time consumption. Through stress tests on multiple sets of mining area data and different grid sizes, setting the checkpoint budget to about one-third of the available memory can save enough slices without overflow, and the overall time consumption is minimal and stable. 35% was selected as the engineering value that balances safety margin and efficiency.

[0069] The number of seismic time steps is determined by the record length and sampling interval; the ratio of the number of seismic time steps to the number of savable slices is calculated and rounded up to obtain the checkpoint interval; the slice time index list and the corresponding wavefield slice array are saved; during the adjoint calculation stage, the forward wavefield is recalculated segment by segment using the saved slices as time boundaries, and adjoint source injection and gradient accumulation are performed simultaneously.

[0070] The update direction vector is constructed using a rotating master-driven approach. Specifically, in one outer iteration, only one or two physics pairs are selected to provide update directions for specified parameter components, while the remaining components are set to zero. The update direction vector is obtained by concatenating the gradient vectors corresponding to the selected physics type. The parameter components corresponding to the unselected physics types are set to zero. The gradient vector is calculated by calling the transpose action interface of the sensitivity operator corresponding to the physics type, taking the normalized residual vector as input, and outputting the parameter space vector as the gradient of the physics type.

[0071] The parameter update step size is calculated using interface consistency-driven methods, and is expressed as follows:

[0072] ;

[0073] in, Indicates the parameter update step size. This indicates updating the direction vector. This indicates the positive action interface corresponding to the current main drive physics. This represents the smallest positive number that prevents the denominator from being zero.

[0074] The parameter update step size calculation can be completed simply by calling the interface. The adaptation matrix can be implemented freely. The parameter field before projection is obtained by subtracting the parameter field of the current iteration from the update direction vector according to the parameter update step size.

[0075] It should be noted that the current iteration parameter field is the parameter field state corresponding to the iteration number; the parameter field before projection is the parameter field state obtained by performing an update once based on the current iteration parameter field according to the update direction vector and parameter update step size, and before constraint projection has been performed; the parameter field of the next iteration is the parameter field state obtained after performing constraint projection on the parameter field before projection.

[0076] Furthermore, the pre-projection parameter field is obtained by updating the current iteration parameter field in parallel on the graphics processor according to voxel linear indices. Specifically, for each voxel linear index, the voxel values ​​of the current iteration parameter field's magnetic susceptibility, conductivity, polarization parameter, and velocity are read, as well as the direction components of the update direction vector at the corresponding four components. The parameter update step size is read, and an update is performed once for each of the four types of parameters. The product of the parameter update step size and the direction component is used as the update amount and subtracted from the corresponding voxel value to obtain the pre-projection magnetic susceptibility voxel value array, the pre-projection conductivity voxel value array, the pre-projection polarization parameter voxel value array, and the pre-projection velocity voxel value array. The four pre-projection voxel value arrays together constitute the pre-projection parameter field. The update direction vector and parameter update step size are written to the versioned metadata iteration log field.

[0077] S3. Generate a trigger set based on the acquisition configuration package, and sequentially perform box-constrained projection, total variation projection and wavelet sparse projection on the parameter field before projection to obtain the sparse parameter field. Perform box-constrained projection based on the sparse parameter field to generate the final parameter field.

[0078] Based on the acquisition configuration package, a set of physical observation types is constructed, the data consistency target norm is calculated, and a trigger set is generated. Based on the pre-projection parameter field, combined with the borehole index array and the box-constrained interval table of borehole properties, the voxel values ​​of the pre-projection parameter field at the borehole index position are truncated in parallel to obtain the box-constrained parameter field.

[0079] Furthermore, a set of physical observation types is constructed. Specifically, the observation type identifier field recorded in the magnetic, electrical, induced polarization, seismic, well logging, and borehole sections is read from the acquisition configuration package. For each section, the string value of the observation type identifier field is read, and it is checked whether an index table for the observation type identifier to the data space vector storage location and an index table for the observation type identifier to the standard deviation vector storage location both exist simultaneously. If both exist, the observation type identifier is added to the set of physical observation types. If the corresponding vector index information is missing, the missing item is written to the data missing summary field of the versioned metadata, and the observation type identifier is skipped.

[0080] It should be noted that the sorting of the physical observation type set adopts the order of the serialized fields of the acquisition configuration package as the primary order and the lexicographical order of the observation type identifier strings as the secondary order. The sorted physical observation type set is written into the versioned metadata iteration log field and participates in the subsequent iteration log hash chain to ensure that the same acquisition configuration package corresponds to the same traversal order.

[0081] Iterate through the set of physical observation types one by one. For each type of physical observation, read the corresponding data space vector storage location from the acquisition configuration package; read the length of the one-dimensional array of vectors from the data space vector; calculate the data consistency target norm of the physical observation type and write the target norm into the versioned metadata iteration log field.

[0082] Specifically, the data consistency target norm for physical observation types is expressed as:

[0083] ;

[0084] in, Indicates the first Data consistency target norm for each type of physical observation Indicates the physical observation type index, Indicates the first The length of the data space vector corresponding to each physical observation type.

[0085] It should be noted that each type of physical observation is normalized according to the noise level, and the consistency target norms of various types of data are comparable in terms of dimensions and scale. Setting a safety margin can cover parallel computing, discretization and interpolation errors, as well as a small number of outlier disturbances, thereby improving iterative stability.

[0086] For the same physical observation type, the normalized residual vector L2 norm is read, and compared with the target norm. If the normalized residual vector L2 norm is greater than the target norm, the physical observation type is added to the trigger set; if the normalized residual vector L2 norm is less than or equal to the target norm, the physical observation type is not added to the trigger set. The trigger set is written to the versioned metadata iteration log field as a list of physical observation type identifier strings, and is sorted in the same way as the target norm list.

[0087] Specifically, the drilling index positions are traversed in parallel on the graphics processor according to the drilling index array. The voxel values ​​of magnetic susceptibility, conductivity, polarization parameters, and velocity of the parameter field before projection are read for each drilling index position. The box constraint interval table entries of the drilling properties corresponding to the drilling index positions are read from the acquisition configuration package, and the corresponding lower and upper bounds are read respectively. If the voxel value is less than the lower bound, the voxel value is replaced with the lower bound. If the voxel value is greater than the upper bound, the voxel value is replaced with the upper bound. If the voxel value is between the lower and upper bounds, the voxel value remains unchanged. The box constraint parameter field is obtained after the parallel truncation is completed.

[0088] Calculate the residual ratio, obtain the total variation projection target scalar, perform dual variable iteration on the box-constrained parameter field to generate the total variation parameter field; calculate the unified noise scale based on the trigger set, calculate the soft thresholding threshold, perform three-dimensional orthogonal wavelet transform on the total variation parameter field and perform soft thresholding processing, and then perform inverse wavelet transform to generate the sparse parameter field; based on the sparse parameter field, perform box-constrained projection through the borehole index array and the box-constrained interval table of borehole properties to generate the final parameter field.

[0089] Furthermore, the total variation projection target scalar is determined by referring to the total variation scalar and the residual ratio. The summation range of the residual ratio is consistent with the summation range of the data consistency target norm, both based on the set of physical observation types.

[0090] Specifically, the residual ratio and the total variation projective target scalar are expressed as:

[0091] ;

[0092] ;

[0093] in, Indicates the residual ratio. Represents a set of physical observation types. Indicates the first Normalized residual vectors for each type of physical observation Indicates the first Normalized residual vector L2 norm for each physical observation type Represents the total variation projected target scalar. This indicates a reference to the total variation scalar.

[0094] Furthermore, total variation projection is performed with the box constraint parameter field as input to obtain the total variation parameter field. Total variation projection is implemented using dual variable iteration, and the stopping condition is the square root of the single-precision floating-point machine precision as the convergence threshold. After each round of dual update, the maximum absolute difference between two adjacent rounds of dual variables is calculated. The dual iteration ends when the maximum absolute difference is less than or equal to the square root of the single-precision floating-point machine precision. The entire total variation projection calculation process is executed on the graphics processor, and the total variation projection target scalar, residual ratio, and dual iteration round number are written to the versioned metadata iteration log field.

[0095] Furthermore, sparse projection is performed on the fully varied parameter field to generate a sparse parameter field. The sparse projection is implemented using a three-dimensional orthogonal wavelet transform and soft thresholding process. The wavelet basis used in the three-dimensional orthogonal wavelet transform is obtained by reading the three-dimensional orthogonal wavelet basis type field of the acquisition configuration package.

[0096] The number of multi-scale decomposition layers of the three-dimensional orthogonal wavelet transform is obtained by calculating the dimension of the three-dimensional voxel grid. Specifically, the grid dimension is read from the versioned metadata, the maximum number of times the grid dimension can be bisected in the three directions is calculated, the minimum of the maximum number of times in the three directions is taken as the number of multi-scale decomposition layers, and the number of multi-scale decomposition layers is written to the iteration log field of the versioned metadata.

[0097] To obtain the soft thresholding threshold, specifically, read the trigger set from the versioned metadata iteration log field, read the corresponding normalized residual vector item by item according to the physical observation type identifier string in the trigger set, take the absolute value of each element of the normalized residual vector for each type and calculate the median to obtain the robust scale of the normalized residual vector; calculate the median of all robust scales in the trigger set to obtain the uniform noise scale, and write the uniform noise scale into the versioned metadata iteration log field.

[0098] The total number of parameter field coefficients is obtained by calculating the number of voxels in the three-dimensional voxel grid and the number of parameter field components. The number of voxels is obtained by multiplying the grid dimensions. The number of parameter field components consists of four components: magnetic susceptibility parameter field, electrical conductivity parameter field, polarization parameter field, and velocity parameter field. The total number of parameter field coefficients is written to the versioned metadata iteration log field.

[0099] Furthermore, by unifying the noise scale and the total number of parameter field coefficients, the soft threshold is obtained and written into the versioned metadata iteration log field.

[0100] Three-dimensional orthogonal wavelet transforms are performed on the magnetic susceptibility, conductivity, polarization, and velocity parameters of the total variation parameter field to obtain wavelet coefficient arrays. Soft thresholding is then performed on each element of the wavelet coefficient arrays. Inverse three-dimensional orthogonal wavelet transforms are then performed on the soft-thresholded wavelet coefficient arrays to generate sparse parameter fields. On the graphics processor, voxel arrays of magnetic susceptibility, conductivity, polarization, and velocity are generated by voxel linear indexing and written to the versioned metadata iteration log field.

[0101] The soft-thresholding threshold and the soft-thresholding wavelet coefficients are expressed as follows:

[0102] ;

[0103] ;

[0104] in, This indicates the soft threshold. Indicates a uniform noise scale. This represents the total number of parameter field coefficients. This represents the wavelet coefficients after soft thresholding. This represents the wavelet coefficients obtained from the three-dimensional orthogonal wavelet transform of the total variation parameter field.

[0105] Furthermore, box-constrained projection is performed again on the sparse parameter field to obtain the parameter field for the next iteration. The borehole index array and the box-constrained interval table of borehole properties used for box-constrained projection are obtained from the acquisition configuration package. The projection execution method is the same as that for obtaining the box-constrained parameter field by parallel truncation. The projected parameter field covers the voxel values ​​at the borehole index positions in the sparse parameter field.

[0106] The next iteration parameter field is written into the parameter field entry of the versioned metadata iteration log field, and the iteration number is incremented within the same log field. This ensures that the current iteration parameter field entry in the log field corresponds to the updated parameter field state at the start of the next outer iteration. After the start of the next outer iteration, the current iteration parameter field is read from the versioned metadata iteration log field. Using the current iteration parameter field as the multiphysics forward modeling input, normalized prediction vectors for various types of data are generated, and the normalized residual vector is recalculated. Based on the normalized residual vector, the transpose action interface of the sensitivity operator is called to construct the update direction vector, and then the forward action interface of the sensitivity operator is called to calculate the parameter update step size. On the graphics processor, the new pre-projection parameter field is generated in parallel by voxel linear index and enters the constrained projection process.

[0107] The outer iteration stopping conditions include the data consistency criterion and the parameter field convergence criterion. The data consistency criterion is achieved by determining whether the trigger set is empty. Specifically, if the trigger set is empty, it is determined that all physical observation types satisfy the data consistency target norm.

[0108] The convergence criterion for the parameter field is as follows: On the graphics processor, voxel-by-voxel differencing is performed on the four components of the parameter field for the next iteration and the current iteration (magnetic susceptibility, conductivity, polarization, and velocity) to obtain four difference voxel arrays. The absolute value of each voxel in the four difference voxel arrays is then taken, and the maximum value is calculated for all voxel positions within each component to obtain the maximum difference of the four components. Finally, the maximum absolute difference is calculated for the four maximum difference values. Further, the absolute value of each voxel in the four components of the current iteration parameter field is taken, and the maximum absolute difference is calculated for all voxel positions within each component. The maximum value of the element position is taken to obtain the maximum absolute value of the four components; then the maximum value of the maximum absolute value of the four components is taken to obtain the maximum absolute value of the current iteration parameter field; the machine precision constant of IEEE-754 single-precision floating-point type is read, where the machine precision constant is a fixed constant between (0,1), which is obtained by reading through compiler constants or runtime floating-point type constants; the square root of the machine precision constant is calculated to obtain the square root of the machine precision; the convergence threshold is the product of the square root of the machine precision and the amplitude reference, where the amplitude reference is the larger value between the maximum absolute value of the current iteration parameter field and 1.

[0109] Furthermore, when the trigger set is empty and the maximum absolute difference is less than or equal to the convergence threshold, the outer iteration is determined to meet the stopping condition, the next iteration parameter field is marked as the final parameter field, and written to the versioned metadata iteration log field; when the trigger set is not empty or the maximum absolute difference is greater than the convergence threshold, the outer iteration is determined not to meet the stopping condition, and the next round of outer iteration process continues to be executed according to the current iteration parameter field recorded in the versioned metadata iteration log field. Specifically, the next round of outer iteration process involves performing multiphysics forward modeling based on the current iteration parameter field, generating normalized prediction vectors for various types of data, calculating normalized residual vectors, calling the transpose action interface of the sensitivity operator to construct update direction vectors, calling the forward action interface of the sensitivity operator to calculate parameter update step size, generating the pre-projection parameter field, and sequentially performing box-constrained projection, total variation projection, and wavelet sparse projection to generate the next iteration parameter field.

[0110] S4. Based on the final parameter field, generate a three-dimensional grade voxel array, calculate the grade segmentation threshold, extract the ore body connected body set, generate the ore body triangular mesh and the ore body existence probability voxel array, and obtain the digital twin base data package.

[0111] Based on the final parameter field, the borehole index array is scanned in parallel to construct a physical property grade training sample table. Kernel regression is used to establish a grade prediction operator and generate a three-dimensional grade voxel array. The grade segmentation threshold is calculated to generate a binary ore body indicator voxel array. Three-dimensional connected component labeling is performed to extract the ore body connected component set.

[0112] Furthermore, the final parameter field includes four one-dimensional arrays: magnetic susceptibility voxel value array, electrical conductivity voxel value array, polarization parameter voxel value array, and velocity voxel value array. The borehole index array, borehole grade array, box constraint interval table of borehole properties, voxel center coordinate array, and grid step size are read from the acquisition configuration package to establish the spatial correspondence between boreholes and voxel grids.

[0113] Based on the borehole index array, the borehole grade array is mapped to a set of borehole support points in a 3D voxel mesh. The voxel center coordinates of the set of borehole support points are obtained by reading the voxel center coordinate array. A physical property grade training sample table is constructed on the set of borehole support points. Specifically, the physical property grade data is recorded one by one by scanning the borehole index array in parallel. The input fields of the physical property grade training sample table are the voxel values ​​of magnetic susceptibility, conductivity, polarization parameter, and velocity of the final parameter field at the borehole index position. The output field of the physical property grade training sample table is the borehole grade value at the same borehole index position. The physical property grade training sample table is written into versioned metadata and the number of sample entries is recorded.

[0114] Based on the physical property grade training sample table, a grade prediction operator is established. Specifically, the grade prediction operator is implemented using kernel regression. The kernel width is obtained by calculating the pairwise distance distribution of the input fields in the physical property grade training sample table. Specifically, the Euclidean distance is calculated pairwise for each input field vector in the physical property grade training sample table, and the median is used as the kernel width. The regression coefficients of the kernel regression are solved using the conjugate gradient iteration method. The iteration convergence criterion is the square root of the single-precision floating-point machine precision. The single-precision floating-point machine precision is obtained by reading the machine precision constant of the IEEE-754 single-precision floating-point type. The kernel width and convergence criterion are written into the versioned metadata iteration log field.

[0115] Data cleaning and scaling were performed on the training sample table. Specifically, sample records with missing markers in the input or grade fields were removed, quantile truncation was performed, the median and absolute median difference were calculated for the four types of input fields, and the input fields were scaled uniformly according to the median and absolute median difference.

[0116] The training set and the validation set are divided. Specifically, the drill hole identifier is used as the grouping key, the drill hole identifier is mapped to an integer, and the remainder is taken when the remainder is 5. The samples corresponding to the drill holes with a remainder of 0 are assigned to the validation set, and the remaining samples are assigned to the training set.

[0117] Construct a kernel matrix based on the kernel function, where the kernel function is expressed as:

[0118] ;

[0119] in, This represents the kernel similarity between the corresponding input vectors of two sample records. This represents the scaled input field vector of the first sample record. This represents the scaled input field vector of the second sample record. Indicates the kernel width.

[0120] The baseline value for the kernel width is the median of the pairwise Euclidean distances between the input field vectors in the training set. Candidate kernel width values ​​are taken as five fixed proportions of the baseline kernel width: 1 / 2, 3 / 4, the baseline value, 4 / 5, and 2 / 3. These candidate kernel width values ​​are written to the versioned metadata iteration log field. A training set kernel matrix is ​​constructed for each candidate kernel width value, and the regression coefficients are solved iteratively using conjugate gradients. The regularization coefficient and the conjugate gradient stopping threshold are both on the order of the square root of single-precision machine precision. The stopping threshold is determined by combining the L2 norm of the training set quality vectors. The regularization coefficient, stopping threshold, maximum number of iterations, and actual number of iterations are written to the versioned metadata iteration log field.

[0121] Calculate the root mean square error on the validation set, select the kernel width corresponding to the minimum error as the kernel width of the grade prediction operator, and write the error value into the versioned metadata iteration log field.

[0122] Parallel grade prediction is performed at each voxel center of the 3D voxel mesh using a grade prediction operator. The grade prediction inputs are the voxel values ​​of magnetic susceptibility, conductivity, polarization parameters, and velocity of the final parameter field at the voxel center. The grade prediction output is the voxel grade value, generating a 3D grade voxel array. The 3D grade voxel array is then written into versioned metadata and a digital twin base data package.

[0123] Calculate the grade segmentation threshold on the borehole grade array. The grade segmentation threshold is obtained by creating a histogram of borehole grade samples and performing inter-class variance maximization segmentation. Write the grade segmentation threshold into versioned metadata.

[0124] It should be noted that, specifically, obtaining the grade segmentation threshold involves reading the borehole grade array from the acquisition configuration package, generating a one-dimensional sample set according to the grade values ​​corresponding to the borehole index array, removing samples with missing markers, and writing the number of missing markers into the versioned metadata; calculating the minimum and maximum values ​​of the sample set as the lower and upper bounds of the range for the grade segmentation threshold search, and writing the lower and upper bounds into the versioned metadata; constructing a fixed-bin histogram based on the range, with the number of bins obtained by rounding up the square root of the number of samples, and the histogram bin boundaries being generated at equal intervals between the lower and upper bounds; calculating the frequency proportion and mean of samples on both sides of each candidate grade segmentation threshold, and selecting the candidate grade segmentation threshold with the largest inter-class variance as the grade segmentation threshold using the inter-class variance maximization criterion; the range of the grade segmentation threshold is between the minimum and maximum values ​​of the sample set, and the grade segmentation threshold is one of the bin boundary points of the histogram.

[0125] Based on a three-dimensional grade voxel array and a grade segmentation threshold, a binary ore body indicator voxel array is generated by comparing the grade value at each voxel location with the grade segmentation threshold.

[0126] Three-dimensional connected component labeling is performed on the ore body indicator voxel array to extract the ore body connected component set. The connected component labeling adopts a six-adjacency topology, which is determined by the voxel linear index and three-dimensional index mapping array. For each connected component in the ore body connected component set, its boundary voxels are extracted, and the corresponding ore body triangular mesh is generated through the isosurface extraction process. The isosurface values ​​are determined by the binary definition of the ore body indicator voxel array. The ore body triangular mesh is written into the digital twin base data package.

[0127] For each orebody in the set of orebody interconnections, extract the boundary voxels of the orebody interconnections to generate an orebody triangular mesh; based on the standard deviation vector, construct a set of perturbation observation vectors, perform a single round of outer layer iteration to obtain an array of voxels representing the orebody existence probability, and combine the three-dimensional grade voxel array and the orebody triangular mesh to obtain a digital twin base data package.

[0128] Furthermore, a random perturbation is applied to the normalized observation vector using the standard deviation vector to construct a perturbed observation vector set. The random number seed is obtained by truncating the low-order bits of the versioned metadata hash value and written into the versioned metadata iteration log field. The number of samples in the perturbed observation vector set is determined by the statistical convergence criterion. Specifically, a single-round outer iteration is performed on each generated perturbed observation vector to generate a corresponding ore body indicator voxel array. The occurrence frequency of the ore body indicator voxel array at each voxel position is calculated cumulatively. After each round of updates, the maximum absolute difference between the occurrence frequency array and the occurrence frequency array of the previous round is calculated. When the maximum absolute difference is less than or equal to the square root of the single-precision floating-point machine precision, the generation of perturbed observation vectors is stopped. The final occurrence frequency array is used as the ore body existence probability voxel array and written into the digital twin base data package.

[0129] The digital twin base data package is serialized in a fixed field order. The digital twin base data package includes, but is not limited to, the three-dimensional voxel grid definition, the final parameter field four types of physical property voxel value array, the three-dimensional grade voxel array, the ore body triangular grid, the ore body existence probability voxel array, versioned metadata, and the versioned metadata iteration log field; the serialized text participates in the versioned metadata hash link.

[0130] In this embodiment, by comparing different parallel parameters and iteration process logs, the effects of the GPU parallel computing framework on improving the speed of large-scale data processing, shortening computation time, reducing storage pressure, improving processing efficiency, and enhancing convergence controllability and inversion stability in the 3D ore body inversion process were verified. The computing platform adopts a combination of a graphics processor and a host processor with CUDA capabilities, and is equipped with CUDA drivers and runtime environment. Timing adopts GPU-side event timing or equivalent timing method, and memory usage adopts memory peak sampling or equivalent sampling method. The running log and versioned metadata iteration log fields are written to solid-state storage.

[0131] Furthermore, the parallel verification section selects the 3D voxel mesh size as the scale axis, and sets the number of elements within a data block and the number of concurrent processing blocks as the parallelism axes. The number of concurrent processing blocks reflects the concurrency intensity of processing multiple data blocks simultaneously, while the number of elements within a data block reflects the granularity of data block partitioning. Each set of parameters is repeatedly executed multiple rounds of outer layer iterations, and the average time of a single outer layer iteration is recorded and summarized in the form of mean and standard deviation. Figure 2 This is a graph showing the average time taken for a single outer layer iteration versus the total number of voxels. The legend shows different curves corresponding to different numbers of concurrent processing blocks. Figure 2The data shows that as the total number of voxels increases, the average time for a single outer iteration increases, reflecting the increased computational load of large-scale cases. Under the same total number of voxels, as the number of concurrent processing blocks increases, the average time for a single outer iteration decreases, and the difference in curves is more obvious in large-scale intervals. This reflects that the GPU parallel computing framework can more fully improve throughput under large-scale data conditions, and has the effect of shortening computation time and improving processing efficiency.

[0132] Furthermore, the convergence controllability and inversion stability verification uses the same acquisition configuration package and the same normalization caliber. During the outer iteration process, the trigger set size, unified noise scale, and soft thresholding threshold are recorded synchronously. The normalized residual vector L2 norm of each observation type is used to form a time series. Figure 3 This is a graph showing the norm of the normalized residual vector versus the outer iteration number for each observation type. In the legend, each curve corresponds to a different observation type identifier string. Figure 3 The normalized residual vector L2 norm of multiple observation types shows an overall decrease with the outer iteration number, and the curve shape is continuous and smooth without significant oscillations or repeated jumps. This indicates that the projection constraint chain reduces the residuals of each observation type in a controllable manner while fusing prior information and observation data. This embodiment also provides a computer device suitable for GPU-based rapid inversion methods of 3D ore bodies, including: a memory and a processor; the memory stores computer-executable instructions, and the processor executes the computer-executable instructions to realize the GPU-based rapid inversion method of 3D ore bodies proposed in the above embodiment.

[0133] The computer device can be a terminal, comprising a processor, memory, communication interface, display screen, and input devices connected via a system bus. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, carrier networks, NFC (Near Field Communication), or other technologies. The display screen can be an LCD screen or an e-ink screen. The input devices can be a touch layer covering the display screen, buttons, a trackball, or a touchpad on the computer device's casing, or an external keyboard, touchpad, or mouse.

[0134] This embodiment also provides a storage medium storing a computer program that, when executed by a processor, implements the GPU-based fast inversion method for 3D ore bodies as proposed in the above embodiments. The storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read Only Memory (EPROM), Programmable Red-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk.

[0135] In summary, this invention improves the speed of large-scale data processing in the 3D ore body inversion process by simultaneously processing multiple data blocks using a GPU parallel computing framework, thereby shortening computation time, reducing storage pressure, and increasing processing efficiency. Furthermore, it enhances convergence controllability and inversion stability by accurately fusing prior information and observation data through an optimization algorithm based on projection constraints and by adaptively determining the soft thresholding threshold based on trigger sets and a unified noise scale.

[0136] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A GPU-based method for fast inversion of 3D ore bodies, characterized in that: include, Import multi-source exploration data, generate observation data in a unified coordinate system, establish a three-dimensional voxel grid, generate a borehole index array and a box-constrained interval table of borehole properties, generate an initial parameter field, and obtain the acquisition configuration package. Based on a 3D voxel mesh and acquisition configuration package, multiphysics forward modeling is performed to generate normalized prediction vectors and normalized residual vectors. The forward action interface and transpose action interface of the sensitivity operator are established, the update direction vector is constructed, and the pre-projection parameter field is generated. Based on the acquisition configuration package, a trigger set is generated. Box-constrained projection, total variation projection, and wavelet sparse projection are sequentially performed on the parameter field before projection to obtain the sparse parameter field. Box-constrained projection is then performed on the sparse parameter field to generate the final parameter field. Based on the final parameter field, a three-dimensional grade voxel array is generated, the grade segmentation threshold is calculated, the ore body connected body set is extracted, the ore body triangular mesh and the ore body existence probability voxel array are generated, and the digital twin base data package is obtained. The generation of the trigger set includes: constructing a set of physical observation types based on the acquisition configuration package, calculating the data consistency target norm, and generating the trigger set; The specific steps for generating the final parameter field are as follows: Based on the pre-projection parameter field, combined with the borehole index array and the box-constrained interval table of borehole properties, the voxel values ​​of the pre-projection parameter field at the borehole index position are truncated in parallel to obtain the box-constrained parameter field. Calculate the residual ratio, obtain the total variation projection target scalar, perform dual variable iteration on the box constraint parameter field, and generate the total variation parameter field; Based on the trigger set calculation of the unified noise scale, the soft thresholding threshold is calculated, the three-dimensional orthogonal wavelet transform of the total variation parameter field is performed and soft thresholding is applied, and then the inverse wavelet transform is performed to generate the sparse parameter field. Based on the sparse parameter field, the final parameter field is generated by performing box-constrained projection through the borehole index array and the box-constrained interval table of borehole properties. The specific steps for obtaining the digital twin base data package are as follows: For each orebody in the set of orebody connected bodies, extract the boundary voxels of the orebody connected bodies to generate a orebody triangular mesh; Based on the standard deviation vector, a set of perturbation observation vectors is constructed, a single round of outer layer iteration is performed, an array of voxel arrays of ore body existence probability is obtained, and combined with the three-dimensional grade voxel array and the ore body triangular mesh, a digital twin base data package is obtained.

2. The GPU-based fast inversion method for 3D ore bodies as described in claim 1, characterized in that: The specific steps for obtaining the data acquisition configuration package are as follows: Import multi-source exploration data, calculate coordinate transformation parameters, obtain observation data under a unified coordinate system, and calculate the standard deviation vector of each type of data. Based on observation data in a unified coordinate system, the horizontal and vertical grid step sizes are calculated to generate a three-dimensional voxel grid. In the three-dimensional voxel mesh, the nearest voxel is searched for each borehole trajectory point to generate a borehole index array. Combined with the standard deviation vector, a box-constrained interval table of borehole properties is generated. Based on the borehole index array and the box-constrained interval table of borehole properties, Gaussian radial basis functions are used to interpolate on the three-dimensional voxel mesh to generate the initial parameter field and obtain the acquisition configuration package.

3. The GPU-based fast inversion method for 3D ore bodies as described in claim 2, characterized in that: The specific steps for generating the normalized prediction vector and the normalized residual vector are as follows: Based on the three-dimensional voxel mesh and the acquisition configuration package, multiphysics forward modeling is performed to generate normalized prediction vectors for various types of data. The normalized observation vector and the normalized prediction vector of various types of data are subtracted element by element to obtain the normalized residual vector.

4. The GPU-based fast inversion method for 3D ore bodies as described in claim 3, characterized in that: The process of generating the pre-projection parameter field includes: establishing a positive action interface and a transpose action interface for the sensitivity operator pair vector; calculating the checkpoint interval; calling the transpose action interface based on the normalized residual vector; constructing the update direction vector using a rotating master-driven approach; calling the positive action interface; calculating the parameter update step size; and performing an update along the update direction vector based on the current iterative parameter field to generate the pre-projection parameter field.

5. The GPU-based fast inversion method for 3D ore bodies as described in claim 4, characterized in that: The specific steps for extracting the connected components of the ore body are as follows: Based on the final parameter field, the borehole index array is scanned in parallel to construct a physical property grade training sample table, and a grade prediction operator is established by kernel regression to generate a three-dimensional grade voxel array. Calculate the grade segmentation threshold, generate a binary ore body indicator voxel array, perform three-dimensional connected component labeling, and extract the ore body connected component set.

6. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that: When the processor executes the computer program, it implements the steps of the GPU-based fast inversion method for three-dimensional ore bodies as described in any one of claims 1 to 5.

7. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by the processor, it implements the steps of the GPU-based fast inversion method for three-dimensional ore bodies as described in any one of claims 1 to 5.

Images