Transient electromagnetic inversion method and system with multi-scale physical confidence constraints

Abstract

The application provides a multi-scale physical confidence constraint transient electromagnetic inversion method and system, and the method comprises the following steps: acquiring multi-channel transient electromagnetic data excited by secondary induction electromagnetic field; extracting attenuation characteristics by using a space-time adaptive decoupling network, determining spatial variation confidence by Monte Carlo inference, and constructing a non-uniform adaptive grid and an initial resistivity model according to the spatial variation confidence; converting the confidence into a spatial weight factor and dynamically embedding the spatial weight factor into a damping factor and a regularization matrix of an LM algorithm to construct a target function containing electromagnetic physical constraints; solving a sensitivity matrix by using a multi-field source accompanying state method, calculating a model update amount, iteratively updating the initial model, and outputting a three-dimensional resistivity imaging graph of a spatial physical boundary of a subsurface target body. The application effectively overcomes signal distortion caused by transient electromagnetic wave diffusion scattering, and solves the technical problems of high dependence of a traditional physical inversion on an artificially set initial model and unclear delineation of a subsurface target body physical boundary.

Description

Technical Field

[0001] This invention belongs to the field of electromagnetic data analysis technology, specifically relating to a transient electromagnetic inversion method and system with multi-scale physical confidence constraints. Background Technology

[0002] Transient electromagnetic methods (TEMs) utilize electromagnetic sensors to receive secondary induced electromagnetic fields generated by artificial excitation in underground media. They are a crucial tool for geophysical imaging of underground structures in fields such as underground mineral exploration, coalfield goaf remediation, and engineering geological surveys. In actual exploration, the resistivity distribution of underground rock strata, faults, and water-bearing voids exhibits strong spatial anisotropy and non-uniformity. As transient electromagnetic waves propagate deeper underground, complex physical scattering and boundary reflections occur, resulting in the multi-channel transient electromagnetic signals acquired by the receiving electromagnetic sensors displaying strong and non-stationary time-varying attenuation characteristics along the time axis. Traditional transient electromagnetic inversion typically relies on manually established uniform initial geological resistivity models and iteratively solving Maxwell's equations to reconstruct the underground resistivity distribution.

[0003] However, due to the unknown spatial distribution of the prior electrical conductivity of non-uniform underground media, if the initial resistivity model deviates significantly from the actual geological structure, the electromagnetic forward modeling response in numerical calculations will suffer from severe physical mismatch with the physical signals measured by sensors. This causes the gradient update direction of the inversion objective function to deviate from reality, making the inversion calculation prone to getting trapped in local extrema and resulting in severe distortion of the physical spatial boundary characterization of underground geological bodies. While data-driven neural network methods can establish fast mappings, their outputs do not conform to the constraints of electromagnetic physical equations. When faced with electromagnetic noise in the field and the physical temperature drift of electromagnetic sensors, they are prone to producing spurious electrical interfaces that do not conform to physical laws. Therefore, overcoming signal distortion caused by transient electromagnetic wave diffusion and scattering in non-uniform and complex geological exploration, and solving the problems of high dependence on manually set initial models and unclear characterization of the physical boundaries of underground targets in traditional physical inversion are currently the main challenges. Summary of the Invention

[0004] This invention provides a transient electromagnetic inversion method and system with multi-scale physical confidence constraints to solve the above-mentioned technical problems.

[0005] In a first aspect, the present invention provides a transient electromagnetic inversion method with multi-scale physical confidence constraints, the method comprising the following steps: Multi-channel transient electromagnetic data of underground exploration targets excited by secondary induced electromagnetic fields are obtained by electromagnetic sensors. A spatiotemporal adaptive decoupling network containing a time-series encoder and a spatial decoder is constructed. The time-series encoder adopts a cascaded structure of multi-scale temporal convolution and bidirectional long short-term memory network, and the spatial decoder adopts a spatial transposed convolution structure with an attention mechanism. Multi-channel transient electromagnetic data is input into a spatiotemporal adaptive decoupling network, and multiple forward inferences are performed in the spatiotemporal adaptive decoupling network based on the Monte Carlo random drop mechanism to obtain multiple sets of predicted resistivity distributions. The spatial variation confidence level, which represents the uncertainty of the geomorphological prediction of the underground medium where the underground exploration target is located, is obtained by quantifying the spatial distribution variance of the multiple sets of predicted resistivity distributions. Based on the spatial variation confidence level, a non-uniform adaptive grid is constructed for the underground medium where the underground exploration target is located. Multiple sets of predicted resistivity distributions are discretized in the non-uniform adaptive grid to construct an initial resistivity model. The spatial variation confidence is transformed into a spatial weight factor, and the spatial weight factor is dynamically embedded into the damping factor and prior regularization matrix of the LM algorithm in the physical iterative inversion process to construct an objective function containing physical constraints. The adjoint state method under multiple field source constraints is used to solve the sensitivity matrix of the initial resistivity model in a non-uniform adaptive grid, and the objective function is solved through the sensitivity matrix to obtain the model update amount; The initial resistivity model is iteratively updated using model update increments until the preset convergence condition is met, and a three-dimensional resistivity image of the spatial physical boundary of the underground exploration target is output.

[0006] Optionally, before inputting the multi-channel transient electromagnetic data into the spatiotemporal adaptive decoupling network, the method further includes a step of calculating the vertical magnetic field components of a pre-defined subsurface medium model using transient electromagnetic forward modeling equations as simulated observation data, and using the simulated observation data to construct a training dataset to train the spatiotemporal adaptive decoupling network. The calculation of the vertical magnetic field components using transient electromagnetic forward modeling equations is described in detail below. The calculation formula is as follows: , in, Indicates a dimensionless time parameter. Represents the characteristic time constant. Represents the permeability of free space. Indicates the radius of the transmission loop. Indicates the transmitting current. Represents the resistivity of the subsurface medium model. This indicates the preset decay time. This represents the preset error function.

[0007] Optionally, the construction and output process of the time series encoder includes: Multiple preset temporal convolutional kernels with different receptive fields are configured to construct a multi-scale temporal feature extraction layer; Local multi-scale temporal features of multi-channel transient electromagnetic data are extracted using a multi-scale temporal feature extraction layer. A bidirectional long short-term memory network is constructed to capture the temporal dependencies of local multi-scale temporal features in both the forward and reverse temporal directions, thereby obtaining a global temporal feature vector. Local multi-scale temporal features are concatenated and fused with global temporal feature vectors to obtain concatenated temporal features, which are then used as the output of the time series encoder.

[0008] Optionally, the construction and output process of the spatial decoder includes: Construct multiple spatially transposed convolutional layers to progressively upsample the cascaded temporal features to restore them to three-dimensional spatial features; A dual attention unit consisting of channel attention and spatial attention mechanisms is embedded between adjacent spatially transposed convolutional layers; The attention feature weights of the three-dimensional spatial features in the channel dimension and spatial dimension are calculated using dual attention units; Attention feature weights are used to weight and adjust the three-dimensional spatial features to focus on the boundary features of the resistivity spatial distribution, and multiple sets of predicted resistivity distributions are output.

[0009] Optionally, the step of performing multiple forward inferences in the spatiotemporal adaptive decoupling network based on the Monte Carlo random dropout mechanism to obtain multiple sets of predicted resistivity distributions, and obtaining the spatial variation confidence level representing the uncertainty of the geomorphological prediction of the subsurface medium where the subsurface exploration target is located by statistically quantifying the spatial distribution variance of multiple sets of predicted resistivity distributions, includes the following steps: Configure Monte Carlo random dropout layers in the hidden layers of the spatiotemporal adaptive decoupling network; During the network prediction phase, the Monte Carlo random drop-out layer is kept active. Multi-channel transient electromagnetic data is repeatedly input into the spatiotemporal adaptive decoupling network to perform forward prediction inference a preset number of times, resulting in multiple sets of predicted resistivity distributions. Calculate the average resistivity and spatial distribution variance of multiple predicted resistivity distributions at each spatial node; The reciprocal of the spatial distribution variance is determined as the spatial variation confidence level, representing the uncertainty of predicting the geomorphology of the underground medium where the underground exploration target is located.

[0010] Optionally, the step of constructing a non-uniform adaptive grid of the subsurface medium where the subsurface exploration target is located based on the spatial variation confidence level, and discretizing multiple sets of predicted resistivity distributions in the non-uniform adaptive grid to construct an initial resistivity model includes the following steps: Obtain a preset background spatial grid containing the initial subdivision accuracy, and match the corresponding spatial variation confidence at the nodes of the background spatial grid; In the background spatial grid, identify prediction blind zones where the confidence level of spatial variation is less than a preset confidence threshold, and increase the initial subdivision accuracy within the prediction blind zones; In the background spatial grid, identify confidence regions with spatial variation confidence not less than the confidence threshold, and reduce the initial subdivision accuracy within the confidence regions to construct a non-uniform adaptive grid for the subsurface medium where the subsurface exploration target is located. Multiple sets of predicted resistivity distributions are interpolated and projected onto a non-uniform adaptive grid for discretization to obtain the initial resistivity model.

[0011] Optionally, the step of converting spatial variation confidence into spatial weighting factors and dynamically embedding spatial weighting factors into the damping factor and prior regularization matrix of the LM algorithm in the physical iterative inversion process to construct an objective function containing physical constraints includes the following steps: The spatial variation confidence level at each grid cell in the initial resistivity model is determined as the spatial weighting factor. Construct a diagonal matrix to constrain the spatial inversion weights using spatial weighting factors; By using a diagonal matrix to perform a product weighting adjustment on a pre-defined prior regularization matrix, the regularization constraint strength of a non-uniform distribution can be obtained. The preset damping factor matrix of the LM algorithm during the iteration process is adjusted by using a diagonal matrix to increase the adaptive damping coefficient at the prediction blind zone, so as to obtain the adjusted damping factor matrix. The regularization constraint strength and the adjusted damping factor matrix are used as physical confidence constraints and embedded in the physical iteration of the LM algorithm to construct the objective function.

[0012] Optionally, the objective function formula is as follows: , in, This represents the value of the objective function. This represents the current model parameters after discretization on a non-uniform adaptive grid. This represents multi-channel transient electromagnetic data. Indicates the current model parameters Forward modeling composite data on a non-uniform adaptive grid Represents the adaptive regularization matrix. This represents the preset regularization parameters.

[0013] Optionally, the model parameter update equation corresponding to the objective function is specifically as follows: , in, Indicates the first The sensitivity matrix corresponding to the next iteration Indicates the first The amount of model parameter updates in each iteration. Indicates the first Model parameters for the next iteration Indicates the first The forward-modeled composite data of the model parameters in the next iteration on a non-uniform adaptive grid. This represents the preset damping control coefficient. This represents the adjusted damping factor matrix.

[0014] In a second aspect, the present invention also provides a transient electromagnetic inversion system with multi-scale physical confidence constraints, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the transient electromagnetic inversion method with multi-scale physical confidence constraints as described in any one of the first aspects.

[0015] The beneficial effects of this invention are: This invention utilizes multi-channel transient electromagnetic data acquired from underground exploration targets using electromagnetic sensors. It employs a spatiotemporal adaptive decoupling network to physically decouple the temporal decay characteristics and spatial resistivity distribution of transient electromagnetic signals, enabling rapid inversion of the electrical characteristics of underground media at different depths from highly decaying transient signals. Based on this, multiple forward inferences are performed using a Monte Carlo mechanism to statistically and quantify the spatial variation confidence level of the uncertainty in predicting underground media properties. This physically and quantitatively characterizes the measurement uncertainty of electrical signals caused by random electromagnetic interference in the field and physical temperature drift of the sensors. Based on this, a non-uniform adaptive grid for the underground media is constructed, dynamically matching the spatial grid density with the severity of abrupt changes in underground electrical properties. Then, the spatial confidence level is transformed into a spatial weighting factor, which is dynamically embedded into the damping factor and regularization matrix of the physical iterative inversion. This ensures that the solution of the objective function is directly constrained by the boundary constraints of the physical equations for the diffusion of underground electromagnetic waves. The adaptive spatial adjustment of the damping factor guarantees that the inversion update direction always follows the physical solution space that conforms to the conservation of electromagnetic field energy. This eliminates the absolute dependence of the inversion process on the initial model set by human experience, effectively avoiding physical boundary ambiguity and spurious electrical distortion caused by initial value deviations. Finally, the adjoint state method under multi-source constraints is used to quickly solve the sensitivity matrix in the non-uniform grid, greatly reducing the physical overhead of the forward modeling of three-dimensional partial differential equations. While ensuring high efficiency, this significantly improves the spatial resolution and geometric boundary positioning accuracy of the three-dimensional resistivity imaging of the spatial physical boundary of the underground exploration target. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating a transient electromagnetic inversion method with multi-scale physical confidence constraints in one embodiment of this application.

[0017] Figure 2 This is a data comparison graph showing the attenuation voltage curves of different methods and actual observation data in one embodiment of this application.

[0018] Figure 3 This is a data comparison diagram of different methods in the iterative calculus of the Levenburg-Mcquart algorithm in one embodiment of this application. Detailed Implementation

[0019] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0020] The terms "first," "second," etc., used in the specification and claims of this application are used to distinguish similar objects and not to describe a specific order or sequence. It should be understood that such use of data can be interchanged where appropriate so that embodiments of this application can be implemented in orders other than those illustrated or described herein, and the objects distinguished by "first," "second," etc., are generally of the same class and the number of objects is not limited; for example, a first object can be one or more. Furthermore, in the specification and claims, "and / or" indicates at least one of the connected objects, and the character " / " generally indicates that the preceding and following objects are in an "or" relationship.

[0021] Figure 1 This is a flowchart illustrating a transient electromagnetic inversion method with multi-scale physical confidence constraints in one embodiment. It should be understood that, although... Figure 1 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified in this document, there is no strict order in which these steps are executed, and they can be performed in other orders. Furthermore, Figure 1 At least some steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily executed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed alternately or in turn with other steps or at least a portion of the sub-steps or stages of other steps. For example Figure 1 As shown, the transient electromagnetic inversion method with multi-scale physical confidence constraints disclosed in this invention specifically includes the following steps: S101. Utilize electromagnetic sensors to acquire multi-channel transient electromagnetic data of underground exploration targets excited by secondary induced electromagnetic fields, and construct a spatiotemporal adaptive decoupling network containing a time-series encoder and a spatial decoder.

[0022] The spatiotemporal adaptive decoupling network consists of a time-series encoder and a spatial decoder. The time-series encoder employs a cascaded structure of multi-scale temporal convolution and a bidirectional long short-term memory network. Multiple pre-defined temporal convolutional kernels with different receptive fields are used to construct a multi-scale temporal feature extraction layer. This layer extracts local multi-scale temporal features from the multi-channel transient electromagnetic data. Next, a bidirectional long short-term memory network is constructed to capture the temporal dependencies of the local multi-scale temporal features in both forward and reverse temporal directions, resulting in a global temporal feature vector. The local multi-scale temporal features are then cascaded and fused with the global temporal feature vector to obtain cascaded temporal features, which serve as the output of the time-series encoder. The spatial decoder employs a spatially transposed convolutional structure with an attention mechanism. Multiple spatially transposed convolutional layers are constructed to progressively upsample the cascaded temporal features to reconstruct three-dimensional spatial features. Dual attention units, consisting of channel attention and spatial attention mechanisms, are embedded between adjacent spatially transposed convolutional layers. Dual attention units are used to calculate the weights of 3D spatial features in both the channel and spatial dimensions. The weights are then used to adjust the features to focus on the resistivity boundary, outputting multiple sets of predicted resistivity distributions.

[0023] The time-series encoder consists of a multi-scale one-dimensional temporal convolutional layer and a bidirectional long short-term memory network layer connected in series. The one-dimensional temporal convolutional layer uses three temporal convolutional kernels of different sizes to extract local features across different time spans. The spatial decoder consists of three spatially transposed convolutional layers and embedded channel-space dual attention units. Before training the network, the vertical magnetic field components of a pre-defined subsurface medium model are calculated using transient electromagnetic forward modeling equations as simulated observation data to construct the training dataset. The calculation of the vertical magnetic field components... The formula is: Among them, the dimensionless time parameter is The characteristic time constant is , The permeability of free space, The radius of the emission loop, For transmitting current, For the resistivity of the underground medium model, For the preset decay time, This is the error function. By changing the resistivity... Design 10,000 different underground medium models, substitute them into the above formula to calculate the corresponding vertical magnetic field components, add noise, pair them with the underground medium models, and send them into the network as a sample set. Use the gradient descent algorithm to adjust the network weight parameters until the training converges.

[0024] The spatiotemporal adaptive decoupling network achieves adaptive feature extraction and uncertainty quantification of spatiotemporal data by deeply integrating with the physical process of transient electromagnetic inversion. In specific applications, the input is set as multi-channel transient electromagnetic data acquired by electromagnetic sensors, and the output is set as multiple sets of predicted resistivity distributions of underground exploration targets. The input multi-channel data is first processed in the time series encoder to extract multi-scale local information in the time dimension through temporal convolutional kernels with different receptive fields. Then, a bidirectional long short-term memory network is used to extract the global attenuation law of electromagnetic signals from two consecutive time directions, thereby transforming the one-dimensional time response into cascaded temporal features. The output cascaded temporal features are then fed into the spatial decoder. During the process of transposed convolution upsampling to restore the spatial dimension, attention weights are applied to the three-dimensional spatial features through dual attention units of channels and space, so that the network focuses on the boundary regions where resistivity changes sharply, improving the quality of boundary reconstruction. Meanwhile, the network performs repeated predictions through an internally configured Monte Carlo random dropout layer, defining the discrete distribution variance of the multiple predictions as the spatial variation confidence of uncertainty. This serves as the input for subsequent non-uniform grid adaptive partitioning and dynamic updates of the damping factor and regularization matrix of the LM algorithm, thus establishing a link between artificial intelligence prediction and physical inversion.

[0025] S102. Input multi-channel transient electromagnetic data into a spatiotemporal adaptive decoupling network, and perform multiple forward inferences in the spatiotemporal adaptive decoupling network based on the Monte Carlo random drop mechanism to obtain multiple sets of predicted resistivity distributions. By statistically quantifying the spatial distribution variance of multiple sets of predicted resistivity distributions, the spatial variation confidence level representing the uncertainty of the geomorphological prediction of the underground medium where the underground exploration target is located is obtained.

[0026] Multi-channel transient electromagnetic data is input into a spatiotemporal adaptive decoupling network. A Monte Carlo random drop-out layer is configured in the hidden layer of the spatiotemporal adaptive decoupling network. During the network prediction phase, the Monte Carlo random drop-out layer is kept active, and the multi-channel transient electromagnetic data is repeatedly input into the spatiotemporal adaptive decoupling network to perform 100 forward prediction inferences, resulting in 100 sets of predicted resistivity distributions. The average resistivity and spatial distribution variance of the 100 sets of predicted resistivity distributions at each spatial node are calculated. The spatial distribution variance reflects the fluctuation of the network during the forward inference process, reflecting the differences in prediction results caused by network structure uncertainties. The reciprocal of the spatial distribution variance is determined as the spatial variation confidence level. Since a larger variance indicates higher prediction uncertainty, a lower reciprocal of the variance, i.e., a lower spatial variation confidence level, indicates a lower reliability in predicting the geomorphology of the subsurface medium where the subsurface exploration target is located. Conversely, a smaller variance corresponds to a higher spatial variation confidence level, representing a more stable and reliable prediction result.

[0027] S103. Construct a non-uniform adaptive grid for the underground medium where the underground exploration target is located based on the spatial variation confidence level, and discretize multiple sets of predicted resistivity distributions in the non-uniform adaptive grid to construct an initial resistivity model.

[0028] The process involves acquiring a pre-defined background spatial grid with initial subdivision accuracy and matching corresponding spatial variation confidence levels at the nodes of the background spatial grid. Prediction blind zones with spatial variation confidence levels below a pre-set confidence threshold are identified within the background spatial grid. Within these blind zones, due to the low network prediction confidence, the initial subdivision accuracy is increased by refining the grid to improve inversion resolution. Confidence regions with spatial variation confidence levels not less than a confidence threshold are identified within the background spatial grid. Within these confidence regions, the initial subdivision accuracy can be appropriately reduced to coarsen the grid, thus reducing subsequent numerical computation. Through this dynamic adjustment of grid subdivision accuracy, a non-uniform adaptive grid for the subsurface medium where the subsurface exploration target is located is constructed. Finally, an interpolation algorithm is used to project multiple sets of predicted resistivity distributions onto the non-uniform adaptive grid for discretization, calculating the initial resistivity values ​​at the nodes of the non-uniform grid, thereby completing the construction of the initial resistivity model.

[0029] S104. The spatial variation confidence is transformed into a spatial weight factor, and the spatial weight factor is dynamically embedded into the damping factor and prior regularization matrix of the LM algorithm in the physical iterative inversion process to construct an objective function containing physical constraints.

[0030] In this approach, the spatial variation confidence level at each grid cell in the initial resistivity model is directly determined as the spatial weighting factor. A diagonal matrix is ​​constructed using this spatial weighting factor to constrain the spatial inversion weights. The diagonal matrix is ​​then used to perform a product-weighted adjustment on the pre-defined prior regularization matrix to obtain a non-uniformly distributed regularization constraint strength. The diagonal matrix is ​​then used to adjust the pre-defined damping factor matrix used in the LM algorithm's iteration process, increasing the adaptive damping coefficient at the prediction blind zone to obtain the adjusted damping factor matrix. The regularization constraint strength and the adjusted damping factor matrix are then used as physical confidence constraints, embedded in the physical iteration of the LM algorithm, to construct the objective function. The objective function expression is: In the formula, This represents the value of the objective function. This represents the current model parameters after discretization on a non-uniform adaptive grid. This represents multi-channel transient electromagnetic data. This represents the forward composite data of the current model parameters on a non-uniform adaptive grid. Represents the adaptive regularization matrix. This represents the preset regularization parameters.

[0031] S105. Solve the sensitivity matrix of the initial resistivity model in the non-uniform adaptive grid using the adjoint state method under multiple field source constraints, and solve the objective function through the sensitivity matrix to obtain the model update amount.

[0032] In the calculation process, the adjoint state method is used to efficiently solve the sensitivity matrix on the reconstructed non-uniform adaptive mesh, avoiding the computational burden of directly solving partial derivatives. The specific solution process is as follows: The forward modeling governing equations for the subsurface medium in the time domain are established on a non-uniform adaptive grid, and their Maxwell equations take the following form: , , In the formula, and These represent the positive electric field and the positive magnetic field, respectively. Represents the permeability of free space. This represents the excited source current density. It represents the electrical conductivity of the underground medium.

[0033] Using inversion observation residuals As a virtual companion source By applying reverse loading at the corresponding observation points, the corresponding reverse time propagation adjoint equation is constructed: , , In the formula, and They represent the accompanying electric field and the accompanying magnetic field, respectively. denoted by , where represents the reverse propagation time, and T represents the observation time of maximum attenuation.

[0034] Finally, the resistivity of the j-th grid cell in the non-uniform adaptive grid is calculated using the spatiotemporal coincidence inner product integral of the forward electric field and the accompanying electric field. Corresponding sensitivity matrix element : In the formula, This represents the volume space of the j-th non-uniform grid cell. This represents the observation data of the i-th channel. By traversing all channels and grid cells, a complete sensitivity matrix is ​​quickly constructed. .

[0035] The model parameter update equation corresponding to the objective function is calculated using the following formula: In the formula, This represents the sensitivity matrix corresponding to the k-th iteration. This represents the amount of model parameter updates in the k-th iteration. This represents the model parameters in the k-th iteration. This represents the forward composite data of the model parameters in the k-th iteration on a non-uniform adaptive grid. This represents the preset damping control coefficient. This represents the adjusted damping factor matrix. By solving the above linear equations, the model update can be directly calculated, ensuring not only the mathematical stability of the physical solution process but also effectively utilizing the previously extracted confidence information.

[0036] S106. Iteratively update the initial resistivity model using the model update amount until the preset convergence condition is met, and output a three-dimensional resistivity image of the spatial physical boundary of the underground exploration target.

[0037] The process involves iteratively updating the initial resistivity model using model update increments until a preset convergence condition is met, outputting a 3D resistivity image of the spatial physical boundary of the subsurface exploration target. During each iteration, the calculated model update increments are accumulated into the model parameters of the current iteration, generating an updated resistivity model. Next, forward modeling is performed using the updated resistivity model to obtain new forward-modeled composite data. The fitting residual between the new forward-modeled composite data and the multi-channel transient electromagnetic data is calculated, and it is determined whether the fitting residual is less than a preset residual threshold of 0.01, or whether the number of iterations has reached a preset maximum of 50 iterations. If the convergence condition is not met, the sensitivity matrix is ​​solved and new model update increments are calculated, repeating the iterative update process. If the convergence condition is met, the iteration stops, and the latest resistivity model at this point is used as the final inversion result. The final resistivity model is then visualized in 3D on a non-uniform adaptive grid, outputting a 3D resistivity image of the spatial physical boundary of the subsurface exploration target. The imaging clearly shows the boundaries and distribution range of underground low-resistivity or high-resistivity bodies, realizing transient electromagnetic inversion under multi-scale physical confidence constraints.

[0038] In one implementation, the construction and output process of the time series encoder includes: Multiple preset temporal convolutional kernels with different receptive fields are configured to construct a multi-scale temporal feature extraction layer; Local multi-scale temporal features of multi-channel transient electromagnetic data are extracted using a multi-scale temporal feature extraction layer. A bidirectional long short-term memory network is constructed to capture the temporal dependencies of local multi-scale temporal features in both the forward and reverse temporal directions, thereby obtaining a global temporal feature vector. Local multi-scale temporal features are concatenated and fused with global temporal feature vectors to obtain concatenated temporal features, which are then used as the output of the time series encoder.

[0039] In this embodiment, three one-dimensional temporal convolutional kernels with different receptive field sizes are configured in parallel. Specifically, the first temporal convolutional kernel has a size of 3 to capture rapidly changing local details in the early decay of transient electromagnetic signals; the second temporal convolutional kernel has a size of 5 to capture the transitional fluctuations in the mid-term decay signal; and the third temporal convolutional kernel has a size of 7 to capture the slowly changing background long-term trend information in the late decay signal. By configuring multi-scale temporal convolutional kernels in parallel, time series can be scanned simultaneously at small, medium, and large receptive fields, thereby establishing a feature detector that can adaptively match different decay stages. Each scale of the one-dimensional temporal convolutional layer is configured with the same number of output channels, for example, 64 channels, and a zero-padding strategy is used in the convolution operation to ensure that the features extracted by different convolutional kernels remain completely consistent in terms of time step and temporal length. By combining the three parallel convolutional structures with complementary receptive fields in a horizontal dimension, a multi-scale temporal feature extraction layer capable of capturing the temporal evolution characteristics of transient electromagnetic signals in a comprehensive and multi-granular manner is finally formed.

[0040] The acquired multi-channel transient electromagnetic data is input into the constructed multi-scale temporal feature extraction layer. During feature extraction, the multi-channel transient electromagnetic data is treated as a two-dimensional matrix, with the time dimension representing the number of time sampling points and the channel dimension corresponding to the number of channels in the multi-channel electromagnetic receiver sensor. Three one-dimensional temporal convolutional kernels with different receptive fields slide along the time dimension, performing sliding inner product operations with the corresponding local input data, and filtering is then performed at different temporal receptive fields. The specific extraction calculation formula is as follows: In the above formula, Indicates the first The scale in the th... Local features at each time step Indicates the first The convolution kernel at each scale in the first... The coefficient at each weight position Indicates the first Multichannel transient electromagnetic data at each time step Indicates the first The kernel size can be set to 3, 5, or 7. Indicates the first Bias at each scale This represents the activation function. After the one-dimensional temporal features at each scale are extracted, a non-linear mapping is performed using the activation function to eliminate background noise. Finally, the features are merged and concatenated along the channel dimension, while maintaining the total time series length unchanged before and after concatenation, thus outputting local multi-scale temporal features containing information from multiple receptive fields.

[0041] After acquiring local multi-scale temporal features, a bidirectional long short-term memory (LSTM) network capable of capturing temporal dependencies across the entire time domain needs to be established. In network construction, the bidirectional LSM network consists of two parallel hierarchical structures: a forward LSM network and a reverse LSM network. The forward LSM network layer progressively inputs local multi-scale temporal features in a forward order from the start to the end of time, accumulating the unidirectional temporal influence of early decay features on late decay features in the hidden states. The reverse LSM network layer, on the other hand, scans in a reverse order from the end to the start of time, accumulating the inverse feedback relationship between late electromagnetic responses and early evolution trends in the hidden layers. At each time step, the hidden state vectors of the two opposite directions are concatenated and merged along their corresponding dimensions. After the bidirectional network has processed the sequence data for all time steps, the hidden state of the last forward time step is extracted and combined with the hidden state of the first reverse time step to obtain the final global temporal feature vector representing the global decay trend of the electromagnetic signal.

[0042] After obtaining the multi-scale temporal features representing local variations and the long-range temporal feature vectors representing global variations, a deep feature concatenation and fusion process is required. To avoid losing the abrupt changes in local fluctuations after long bidirectional loop computation, a parallel channel concatenation and fusion mechanism is established. Specifically, the global temporal feature vector is first copied and expanded along the time axis through a broadcast mechanism, ensuring that its temporal dimension is identical to the sequence length of the local multi-scale temporal features. Subsequently, at each identical time sampling step, the local multi-scale temporal feature vector of the current time step is concatenated and combined with the expanded global temporal feature vector along the channel dimension. By combining local fine-tuning with the global framework in a concatenation manner, the feature matrix retains the multi-scale high-frequency details of local fluctuations at different sampling points while deeply integrating the global overall evolution state of electromagnetic decay. Finally, the resulting concatenated temporal features are passed as the output of the time series encoder to the spatial decoder.

[0043] In one implementation, the construction and output process of the spatial decoder includes: Construct multiple spatially transposed convolutional layers to progressively upsample the cascaded temporal features to restore them to three-dimensional spatial features; A dual attention unit consisting of channel attention and spatial attention mechanisms is embedded between adjacent spatially transposed convolutional layers; The attention feature weights of the three-dimensional spatial features in the channel dimension and spatial dimension are calculated using dual attention units; Attention feature weights are used to weight and adjust the three-dimensional spatial features to focus on the boundary features of the resistivity spatial distribution, and multiple sets of predicted resistivity distributions are output.

[0044] In this embodiment, cascaded temporal features from the time-series encoder are received and transformed into a low-resolution initial 3D spatial tensor through a one-dimensional to three-dimensional reshaping transformation. To progressively improve the resolution of the 3D space, a spatial upsampling network consisting of four cascaded 3D spatial transposed convolutional layers is constructed. Each spatial transposed convolutional layer is configured with specific kernel size, stride, and padding parameters. For example, the kernel size of the first spatial transposed convolutional layer is set to 4x4x4, the stride to 2, and the padding to 1, which is used to simultaneously magnify the low-resolution initial 3D spatial tensor by a factor of 2 in all three spatial dimensions. The subsequent three spatial transposed convolutional layers use similar parameter settings and perform fractional-stride interpolation and convolution operations on the input 3D feature map through deconvolution calculations. During the progressive upsampling process, the number of feature channels gradually decreases from high dimensions, for example, from 256 channels to 128 channels, 64 channels, and 32 channels, while the spatial grid size increases exponentially. This multi-level step-by-step upsampling method can smoothly reconstruct the three-dimensional spatial topology, effectively alleviating the information discontinuity caused by direct high-magnification upsampling, thereby restoring the cascaded temporal features into three-dimensional spatial features containing underground spatial distribution features.

[0045] Between two adjacent spatially transposed convolutional layers, dual attention units are sequentially embedded to filter and enhance multidimensional electrical features during feature propagation. Each dual attention unit consists of a channel attention submodule and a spatial attention submodule connected in series. The channel attention submodule establishes a correlation model between different feature channels. Since different channels represent the response characteristics of electromagnetic responses at different electrical media and depths, the channel attention mechanism can automatically learn and identify the feature channels most sensitive to subsurface resistivity boundaries, suppressing background noise channels. The spatial attention submodule follows the channel attention submodule, focusing on finding regions of dramatic resistivity changes in three-dimensional space, i.e., the spatial locations of anomaly boundaries and stratigraphic boundaries. By cascading the channel and spatial attention mechanisms, dual adaptive feature reshaping of both channel and spatial dimensions can be achieved as multidimensional spatial features flow through adjacent spatially transposed convolutional layers. This embedding of dual attention units effectively prevents blurring of spatial geometric details during upsampling and enhances the model's resolution in recognizing subsurface electrical abrupt changes.

[0046] When calculating attention weights in the channel and spatial dimensions using a dual attention unit, the intermediate 3D spatial features are first input into the channel attention submodule. Global average pooling and global max pooling are performed in each of the three spatial dimensions to compress the 3D spatial features into two 1D channel description vectors. These two channel description vectors are then fed into a shared multilayer perceptron network for nonlinear transformation. The outputs are summed and processed through an activation function to obtain the channel attention feature weights. The formula for calculating the channel attention feature weights is as follows: In the above formula, This represents the channel attention feature weight vector. This represents the intermediate three-dimensional spatial features of the input. This indicates operations involving global average pooling and multilayer perceptron mapping. This indicates operations involving global max pooling and multilayer perceptron mapping. This represents a preset nonlinear activation function. Next, the features weighted by the channel attention feature weights are subjected to channel average pooling and channel max pooling along the channel axis to obtain two two-dimensional spatial feature maps. These maps are then concatenated and fused using a 3x3x3 three-dimensional convolutional layer. Finally, the spatial attention feature weights are obtained through activation function processing, thereby determining the importance weights of each spatial node in resistivity inversion in three-dimensional space.

[0047] After obtaining the channel attention feature weights and spatial attention feature weights, a weighted adjustment is performed on the 3D spatial features using a product of channel-wise and spatial location-wise operations. First, the channel attention feature weight vector is expanded in dimension and then multiplied element-wise with the input intermediate 3D spatial features, thereby amplifying or suppressing the intensity of each feature channel to obtain channel-reshaped features. Next, the spatial attention feature weight tensor is multiplied element-wise with the channel-reshaped features at their corresponding spatial coordinates to scale the features across different spatial regions. The final weighted feature calculation formula is as follows: In the above formula, This represents the three-dimensional spatial features after dual attention weighting adjustment. Represents the spatial attention feature weights. This represents element-wise multiplication. Through this continuous weighted product, the network automatically filters out irrelevant information such as formation background noise, focusing feature extraction on the boundaries and abrupt changes in the spatial distribution of resistivity. The features output from the last dual attention unit are input into the last spatially transposed convolutional layer, compressing the number of feature channels to 1. After numerical activation mapping, the final output consists of multiple sets of predicted resistivity distributions representing the spatial distribution of resistivity in the subsurface three-dimensional medium.

[0048] In one implementation, multiple forward inferences are performed in a spatiotemporally adaptive decoupling network based on the Monte Carlo random dropout mechanism to obtain multiple sets of predicted resistivity distributions. The spatial variation confidence level, representing the uncertainty of the geomorphological prediction of the subsurface medium where the subsurface exploration target is located, is obtained by quantifying the spatial distribution variance of the multiple sets of predicted resistivity distributions, including the following steps: Configure Monte Carlo random dropout layers in the hidden layers of the spatiotemporal adaptive decoupling network; During the network prediction phase, the Monte Carlo random drop-out layer is kept active. Multi-channel transient electromagnetic data is repeatedly input into the spatiotemporal adaptive decoupling network to perform forward prediction inference a preset number of times, resulting in multiple sets of predicted resistivity distributions. Calculate the average resistivity and spatial distribution variance of multiple predicted resistivity distributions at each spatial node; The reciprocal of the spatial distribution variance is determined as the spatial variation confidence level, representing the uncertainty of predicting the geomorphology of the underground medium where the underground exploration target is located.

[0049] In this embodiment, random dropout operations with specific dropout probabilities are embedded after the bidirectional long short-term memory network layer in the time-series encoder and between adjacent 3D spatial transposed convolutional layers and dual attention units in the spatial decoder. The dropout probability is set to 0.2, meaning that during forward computation, the activation values ​​of neurons in each hidden layer have a 20% probability of being randomly reset to zero, resulting in different network substructures for each computation. Let the feature tensor of the hidden layer be... The output features after processing by the Monte Carlo dropout layer The calculation formula is as follows: In the above formula, This represents the feature tensor after the discarding process. This represents the input hidden layer feature tensor. This represents a random masking factor that follows a Bernoulli distribution with a parameter of 0.8, meaning it takes a value of 0 with a probability of 0.2 and a value of 1 with a probability of 0.8 at each node dimension of the feature. This configuration alters the spatial distribution of the hidden layer feature tensors. By keeping the discarded layers active, the spatiotemporally adaptive decoupling network is transformed into a Bayesian network with probability distribution estimation capabilities, laying a structural foundation for subsequent quantification of uncertainties in subsurface geomorphological predictions.

[0050] When the spatiotemporal adaptive decoupling network enters the prediction and inference stage based on actual exploration data, instead of employing the conventional strategy of disabling random dropout layers, all Monte Carlo random dropout layers in the network are forcibly set to an active state. Pre-processed multi-channel transient electromagnetic data is used as the input signal. While keeping the network parameters constant, the input signal is repeatedly and continuously fed into the spatiotemporal adaptive decoupling network, independently executing 100 forward prediction and inference calculations. Because the Monte Carlo random dropout layers are always active, the set of neurons that are disabled within the network is different each time a forward prediction and inference occurs. This process is equivalent to using 100 network sub-models with slight structural differences in 100 calculations. Let the... The formula for predicting the resistivity distribution is as follows: In the above formula, Indicates the first The predicted resistivity distribution obtained from the forward prediction inference. This represents the input multi-channel transient electromagnetic data. This represents a spatiotemporal adaptive decoupling network in a discarded activation state. Indicates the first The random mask set generated during the forward inference process. Through repeated prediction, 100 sets of three-dimensional predicted resistivity distributions with the same spatial topology but slight differences in the values ​​of each node are finally calculated and output, thus transforming the deterministic single electrical distribution prediction result into a set of samples containing spatial statistical characteristics.

[0051] After obtaining 100 sets of predicted resistivity distributions, statistical indices need to be calculated for each node location in the 3D spatial grid. Each spatial node in the 3D grid is traversed, and resistivity values ​​corresponding to the same node location in 100 sets of models are extracted, forming a statistical sample array containing 100 resistivity values. The average resistivity at each spatial node is calculated, and the obtained average resistivity is used as the expected predicted value for the corresponding node. Subsequently, for each spatial node, the spatial distribution variance of these 100 values ​​at the corresponding node location is further calculated to quantitatively evaluate the dispersion of the deep learning network's predictions at different underground locations. The formula for calculating the spatial distribution variance is as follows: In the above formula, Indicates the first Spatial distribution variance at each spatial node Indicates the first The first group of predicted resistivity distributions The resistivity values ​​at each spatial node. Indicates the first Average resistivity at each spatial node This represents the total number of forward prediction inferences, with a value of 100. When the spatial distribution variance at a certain spatial node is small, it indicates that the electrical prediction results of the corresponding node are highly consistent in the forward prediction of the neural network with multiple structural changes, and the prediction stability is high. Conversely, if the variance is large, it indicates that the prediction results of the network at the corresponding spatial node fluctuate significantly. By calculating the variance of all node features, the spatial variance feature distribution of the entire set of data is directly obtained.

[0052] After calculating the spatial distribution variance of each node in three-dimensional space, in order to establish an evaluation index that can directly guide subsequent physical inversion and has clear physical meaning, it is necessary to convert the spatial distribution variance into spatial variation confidence. This is because the spatial distribution variance directly reflects the dispersion of the network's prediction results at specific spatial nodes and the uncertainty of geomorphic predictions; a larger variance indicates higher prediction uncertainty. To ensure a positive correlation between confidence and numerical value, the inverse of the spatial distribution variance at each spatial node is calculated, thereby obtaining the spatial variation confidence of each node location. The formula for calculating the spatial variation confidence is as follows: In the above formula, Indicates the first Confidence of spatial variation at each spatial node Indicates the first Spatial distribution variance at each spatial node This represents a preset minimum constant to prevent the denominator from being zero. Through the above reciprocal operation, at spatial node locations with large prediction fluctuations and high variance, the corresponding reciprocal value, i.e., the spatial variation confidence, will significantly decrease, accurately outlining the blind zone with high prediction uncertainty. Conversely, at spatial node locations with variance close to zero, the spatial variation confidence will remain high. These characteristics provide a reliable spatial adaptive prior basis for targeted differentiated adaptive grid partitioning and regularized matrix weighting of different regions during subsequent physical inversion.

[0053] In one implementation, a non-uniform adaptive grid for the subsurface medium where the subsurface exploration target is located is constructed based on the spatial variability confidence level. Multiple sets of predicted resistivity distributions are discretized within the non-uniform adaptive grid to construct an initial resistivity model, comprising the following steps: Obtain a preset background spatial grid containing the initial subdivision accuracy, and match the corresponding spatial variation confidence at the nodes of the background spatial grid; In the background spatial grid, identify prediction blind zones where the confidence level of spatial variation is less than a preset confidence threshold, and increase the initial subdivision accuracy within the prediction blind zones; In the background spatial grid, identify confidence regions with spatial variation confidence not less than the confidence threshold, and reduce the initial subdivision accuracy within the confidence regions to construct a non-uniform adaptive grid for the subsurface medium where the subsurface exploration target is located. Multiple sets of predicted resistivity distributions are interpolated and projected onto a non-uniform adaptive grid for discretization to obtain the initial resistivity model.

[0054] In this embodiment, a pre-set background spatial grid with regular cubic subdivisions is obtained. The background spatial grid is uniformly distributed in three-dimensional space, and the subdivision interval in each direction is set to 5 meters, which is the initial subdivision accuracy. Since the spatial node coordinates output by the spatiotemporal adaptive decoupling network completely overlap with the background spatial grid in spatial range, a spatial mapping relationship needs to be established between the two. The specific matching process is as follows: for each node in the background spatial grid, the nearest network output node is found through three-dimensional spatial coordinate retrieval. The spatial variation confidence score calculated at the retrieved network output node is directly assigned to the corresponding background spatial grid node. When the positions of the background spatial grid node and the network output node completely overlap, the corresponding spatial variation confidence score is directly copied to the background spatial grid node. If there is a slight deviation between the two coordinates, a distance-weighted average calculation is performed using the values ​​of neighboring nodes, thereby matching the corresponding spatial variation confidence score value at each node of the entire background spatial grid. This operation, which spatially aligns the uncertainty information of deep learning with the physical grid nodes, can seamlessly import statistical features into the physical grid system, achieving a unified matching of the spatiotemporal adaptive decoupling network output data and the background grid in terms of geometric topology.

[0055] After obtaining the background spatial grid with matched spatial variation confidence, a preset confidence threshold is set. The confidence threshold is typically set to 0.5. The background spatial grid is traversed node by node, and the spatial variation confidence at each grid node is compared with 0.5. If the spatial variation confidence at a grid node is less than 0.5, the corresponding node is determined to be in a region of high prediction uncertainty, and the corresponding spatial range is identified and marked as a prediction blind zone. Prediction blind zones often correspond to areas with abrupt changes in formation electrical properties, greater depth, or lower signal-to-noise ratio. To improve the resolution of electrical boundaries by refining the grid in subsequent physical iteration inversion, a grid refinement operation is performed within the prediction blind zone. Specifically, an octree partitioning algorithm is used to cut each initial grid cell with a side length of 5 meters in the prediction blind zone along the midpoint of the three spatial axes, uniformly dividing it into eight sub-grid cells, reducing the side length of the sub-grids to 2.5 meters.

[0056] During the traversal of the background spatial grid, if the spatial variation confidence value at a certain grid node is found to be no less than the confidence threshold of 0.5, the corresponding node is determined to belong to a region with extremely high reliability in the deep learning model's prediction results, and the corresponding spatial range is identified and marked as a confidence region. Since the spatiotemporal adaptive decoupling network has already provided a highly reliable prediction of the resistivity distribution in the confidence region, there is no need to use an overly fine grid in the subsequent physical inversion to prevent unnecessary computational overhead and multiple solutions. Therefore, a grid coarsening operation is performed within the confidence region. Specifically, a merging algorithm is used to merge adjacent grid cells with a side length of 5 meters, reconstructing them into large grid cells with a side length of 10 meters, thereby reducing the initial subdivision accuracy within the confidence region. By adjusting the grid size in different directions for the prediction blind zone and the confidence region, a non-uniform adaptive grid with a dynamically changing grid size depending on the confidence level is finally constructed throughout the entire underground exploration space. Non-uniform adaptive meshes maintain a dense mesh topology in uncertain regions and a sparse mesh topology in reliable regions, thus ensuring computational efficiency while enabling adaptive resolution allocation for subsequent physical inversion of different terrain prediction regions.

[0057] After successfully reconstructing the non-uniform adaptive mesh, the 100 sets of predicted 3D resistivity distributions calculated using the Monte Carlo random dropout mechanism in the previous steps need to be mapped onto the non-uniform adaptive mesh to construct the physical parameters for subsequent iterations. To maintain the continuity of physical quantities at the mesh nodes, a 3D linear interpolation method is used to project the 100 sets of predicted resistivity distributions onto the newly constructed non-uniform adaptive mesh. The specific resampling calculation formula is as follows: In the above formula, This indicates that the interpolation is mapped to the first non-uniform adaptive grid. The first node at the nth node Group predicted resistivity values, Represents the first in the background space grid The corresponding 8 nodes around the first node Group predicted resistivity values, This represents the corresponding 3D linear interpolation weight coefficients. The above 3D linear interpolation operation is performed at each node location in the non-uniform adaptive mesh. Finally, the average of the 100 predicted resistivity values ​​interpolated on the non-uniform adaptive mesh for each node is calculated. The average value at each node in the non-uniform adaptive mesh is used as the physical model parameter for the corresponding node, thereby constructing an initial resistivity model that includes the prior resistivity distribution.

[0058] In one implementation, the spatial variation confidence is transformed into a spatial weighting factor, and the spatial weighting factor is dynamically embedded into the damping factor and prior regularization matrix of the LM algorithm in the physical iterative inversion process to construct an objective function containing physical constraints, including the following steps: The spatial variation confidence level at each grid cell in the initial resistivity model is determined as the spatial weighting factor. Construct a diagonal matrix to constrain the spatial inversion weights using spatial weighting factors; By using a diagonal matrix to perform a product weighting adjustment on a pre-defined prior regularization matrix, the regularization constraint strength of a non-uniform distribution can be obtained. The preset damping factor matrix of the LM algorithm during the iteration process is adjusted by using a diagonal matrix to increase the adaptive damping coefficient at the prediction blind zone, so as to obtain the adjusted damping factor matrix. The regularization constraint strength and the adjusted damping factor matrix are used as physical confidence constraints and embedded in the physical iteration of the LM algorithm to construct the objective function.

[0059] In this embodiment, the initial resistivity model in the reconstructed non-uniform adaptive grid is used as the basis, and the spatial variation confidence level associated with each independent grid cell is extracted. Since the non-uniform adaptive grid has different partitioning sizes in the prediction blind zone and the confidence region, each grid cell directly corresponds to a specific spatial variation confidence level value. To directly introduce the uncertainty measure derived from the variance of the deep learning network output into the physical inversion, the spatial variation confidence level at each grid cell is directly determined as the corresponding spatial weight factor. The physical meaning of the spatial weight factor is to quantitatively characterize the degree of confidence the physical inversion algorithm has in the electrical parameters at each spatial location in the initial resistivity model. When the spatial variation confidence level of a certain grid cell is high, the corresponding spatial weight factor value is close to 1, indicating a high degree of confidence in the initial resistivity prediction result for the corresponding region; when a certain grid cell is in the prediction blind zone, the corresponding spatial variation confidence level is low, and the corresponding spatial weight factor value approaches 0, indicating an extremely low degree of confidence in the initial resistivity prediction result for the corresponding region.

[0060] After determining the spatial weight factors for all grid cells in the non-uniform adaptive mesh, a diagonal matrix is ​​constructed using these factors to constrain the spatial inversion weights. The dimension of the diagonal matrix is ​​kept exactly the same as the total number of model parameters used in the inversion within the non-uniform adaptive mesh. In practice, the diagonal matrix is ​​constructed as follows: the spatial weight factors corresponding to each grid cell are arranged sequentially along the main diagonal of the diagonal matrix according to their topological index order within the non-uniform adaptive mesh, while all off-diagonal elements are filled with zeros. Through this design, each row of the diagonal matrix corresponds to a model parameter within the non-uniform adaptive mesh. The mathematical expression for the diagonal matrix is ​​as follows: In the above formula, This represents the diagonal matrix used for the weights in the constraint space inversion. In a non-uniform adaptive mesh, the first... Spatial weight factor corresponding to each grid cell This represents the total number of grid cells in a non-uniform adaptive grid. Using a diagonal matrix, the uncertainties in terrain prediction in three-dimensional space can be preserved, and these uncertainties can be directly applied to the physical inversion calculation through matrix multiplication.

[0061] After obtaining the diagonal matrix used to constrain the spatial inversion weights, a pre-defined prior regularization matrix is ​​adjusted using the diagonal matrix through a product weighting process. The prior regularization matrix is ​​used to constrain the spatial roughness of the inversion model. To achieve differentiated spatial constraints, the constructed diagonal matrix is ​​multiplied on the left and right by the pre-defined prior regularization matrix, respectively, so that the regularization strength maps to the spatial weight factors. The weighting adjustment formula for the adaptive regularization matrix is ​​as follows: In the above formula, Represents the adaptive regularization matrix. This represents the diagonal matrix used for the weights in the constraint space inversion. This represents the preset prior regularization matrix. Through algebraic operations of the formula, in the confidence region where the network prediction has high confidence, the spatial weight factor is larger, and the corresponding regularization constraint becomes stronger, forcing the inversion result to be smooth and closely match the predicted value. In the prediction blind region with high uncertainty, the spatial weight factor approaches 0, and the corresponding roughness constraint is extremely small, thereby releasing the physical constraint restrictions and allowing the physical inversion to be corrected on the boundary in a freer physical solution space, eliminating the ambiguity that deep learning may bring and improving imaging accuracy.

[0062] In physical iterative inversion, the Levenberg-McQuart algorithm controls the step size and stability of each iteration through a damping factor. To achieve adaptive control of the physical iteration process, the pre-set damping factor matrix of the Levenberg-McQuart algorithm is dynamically adjusted during the iteration process using a constructed diagonal matrix. In specific implementation, since the spatial weight factor values ​​at the prediction blind zone approach 0, a mapping relationship inversely correlated with the diagonal matrix is ​​constructed to significantly increase the adaptive damping coefficient at the prediction blind zone. This limits the update step size of the blind zone in the early stages of iteration, avoiding numerical divergence and excessive oscillations during the inversion update process. The formula for calculating the adjusted damping factor matrix is ​​as follows: In the above formula, This represents the adjusted damping factor matrix. Represents the identity matrix. This represents the diagonal matrix used for the weights in the constraint space inversion. This represents the preset stability constant. By adjusting the formula, in the prediction blind zone where the spatial weight factor is very small, the adaptive damping coefficient is increased, making the nonlinear solution update more robust in the uncertain region and ensuring the smooth convergence of the entire physical inversion numerical solution process.

[0063] The non-uniform distribution regularization constraint strength represented by the adjusted adaptive regularization matrix, along with the adjusted damping factor matrix, are used as physical confidence constraints and dynamically embedded into the physical iterative framework of the Levenberg-McQuart algorithm to construct a physical inversion objective function containing physical constraints. The objective function includes a fitting deviation term between the measured multi-channel transient electromagnetic data and the forward-modeled composite data, as well as a model prior constraint term incorporating the adaptive regularization matrix. By introducing the adaptive regularization matrix during the construction of the objective function, differentiated prior penalties are adaptively applied to different confidence regions during the numerical minimization process.

[0064] In one embodiment, the following is a complete structural description of the spatiotemporal adaptive decoupling network in this application: The spatiotemporal adaptive decoupling network is structurally composed of a cascaded time-series encoder for extracting temporal features and a spatial decoder for reconstructing spatial electrical features. The input multi-channel transient electromagnetic data has a two-dimensional matrix structure with a size equal to the number of time sampling points multiplied by the number of receiving channels. The input of the time-series encoder is directly connected to the multi-channel one-dimensional time-series signal. The multi-scale temporal feature extraction layer serves as the first layer of the encoder, with three one-dimensional temporal convolutional layers of different receptive field sizes connected in parallel. The first one-dimensional temporal convolutional layer has a kernel size of 3, a stride of 1, and padding of 1; the second one-dimensional temporal convolutional layer has a kernel size of 5, a stride of 1, and padding of 2; and the third one-dimensional temporal convolutional layer has a kernel size of 7, a stride of 1, and padding of 3. The number of output channels of the three parallel one-dimensional temporal convolutional layers is set to 64. After performing sliding convolution operation along the time dimension, the output features of the three convolutional layers are merged along the channel dimension by the splicing operator to generate local multi-scale temporal features with 192 channels.

[0065] The output of the concatenation operator is connected to the input of a bidirectional long short-term memory (LSTM) network. The LSM network consists of two layers, with each layer's hidden state dimension set to 64 in both the forward and backward directions. The forward and backward hidden states are concatenated at each time step to obtain a 128-dimensional feature vector. The output of the LSM network is connected to a global temporal feature extraction layer, which extracts features from the forward endpoint and the backward endpoint to form a 128-dimensional global temporal feature vector. This global temporal feature vector is expanded in temporal dimension through a broadcast operation, and then channel-wise concatenated with local multi-scale temporal features to output a 320-dimensional cascaded temporal feature. The output of the cascaded temporal feature is connected to a dimension reshaping layer, which transforms the low-resolution temporal feature matrix into a low-resolution 3D spatial feature tensor of size 4x4x4 with 256 channels through a one-dimensional to three-dimensional transformation.

[0066] The spatial decoder consists of four sequentially connected 3D spatial transposed convolutional layers and three dual attention units sandwiched between adjacent transposed convolutional layers. The first 3D spatial transposed convolutional layer receives the reshaped 3D spatial feature tensor as input. Its kernel size is 4x4x4, stride is 2, padding is 1, and output channels are 128. The output of the first spatial transposed convolutional layer is connected to the first dual attention unit. The channel attention module within the first dual attention unit contains parallel 3D global average pooling and 3D global max pooling layers, which compress the features into one-dimensional channel vectors and connect them to a shared multilayer perceptron with an input dimension of 128, a hidden layer dimension of 8, and an output dimension of 128. The summed outputs of the shared multilayer perceptron are mapped to channel attention feature weights through a non-linear activation function. The spatial attention module contains parallel channel average pooling and channel max pooling layers. The concatenated 2D features are connected to a 3D convolutional layer with a kernel size of 3x3x3, stride of 1, and padding of 1, outputting spatial attention feature weights.

[0067] The input of the second spatially transposed convolutional layer is connected to the output of the first dual attention unit. The parameters of the second spatially transposed convolutional layer are the same as the first layer, but the number of output channels is reduced to 64. Its output is connected to the second dual attention unit. The input of the third spatially transposed convolutional layer is connected to the output of the second dual attention unit, and the number of output channels is reduced to 32. Its output is connected to the third dual attention unit. The input of the fourth spatially transposed convolutional layer is connected to the output of the third dual attention unit, and the number of output channels is set to 1. This allows the spatial decoder to output a three-dimensional matrix containing multiple sets of predicted resistivity distributions.

[0068] The following is the complete training process of the spatiotemporal adaptive decoupling network in this application: Before using the spatiotemporal adaptive decoupling network for subsurface exploration target inversion, complete network parameter training is required. First, a training dataset is constructed, consisting of a large amount of simulated observation data generated through numerical simulation and corresponding subsurface resistivity distribution labels. In practice, 10,000 three-dimensional subsurface medium models containing boundaries of different electrical anomalies and resistivity values ​​of different anomalies are randomly generated. The vertical magnetic field components of the corresponding subsurface medium models are calculated using transient electromagnetic forward modeling equations, serving as simulated observation data. By changing the resistivity... The numerical values ​​were used to generate corresponding multi-channel transient electromagnetic voltage decay curves using the vertical magnetic field component calculation formula. The generated 10,000 sets of vertical magnetic field component data were used as training input samples, and the corresponding three-dimensional underground resistivity model was used as the training label, with the training and validation sets divided in an 8:2 ratio. Next, all trainable weight parameters of the network were initialized using the Kemming normal distribution initialization method, and the bias parameters were initialized to zero. The training hyperparameters were set, with a batch size of 32 and a maximum training epoch of 200 epochs. The Adam optimizer was selected as the algorithm for updating network parameters, with an initial learning rate set to 0.001.

[0069] To ensure smooth convergence of network weights in the later stages of training, a cosine annealing strategy is employed to gradually decay the learning rate to 1 x 10^-6 over 200 epochs. A target loss function is defined to guide network training; this function is a weighted combination of the mean squared error loss function and the 3D structural similarity loss function, thus simultaneously constraining the approximation of resistivity values ​​and the reshaping of spatial geometric boundaries. In each training iteration, a batch of 32 multi-channel transient electromagnetic data is input into the spatiotemporal adaptive decoupling network for forward propagation calculation, yielding the predicted 3D resistivity distribution at the output. The target loss function value between the predicted distribution and the ground truth label is calculated, and the gradient of the target loss function with respect to all convolutional kernel weights, gating weights, and biases in the temporal encoder, spatial decoder, and dual attention unit is calculated using backpropagation along the network structure. The optimizer calculates the value based on the gradient and updates all weight parameters of the network using the gradient. After each training epoch, the trend of the loss function's value is evaluated on the validation set. The network is considered complete when the loss function value on the validation set no longer decreases for 15 consecutive epochs, or when the maximum training steps of 200 epochs are reached. The set of weight parameters with the minimum loss on the validation set is saved as the final structural parameters of the network, completing the full training of the spatiotemporally adaptive decoupled network.

[0070] The following is a detailed implementation process of the solution proposed in this application in this embodiment: Taking the detection of a three-layered geological model with a high-resistivity interlayer at a depth of 200 to 300 meters as an example, pulsed electromagnetic fields are transmitted from the ground to the underground using a transmitter loop. Multi-channel transient electromagnetic attenuation voltage data with a time span of 10 to 40 milliseconds are collected using electromagnetic sensors. At the initial observation stage, i.e., 10 milliseconds, the collected attenuation voltage data is approximately 10^4 mV / m. 2 As time progressed, in the late observation period of 40 milliseconds, the voltage value decayed to 10^-2 mV / m. 2The following is a continuously decaying voltage curve reflecting the electrical distribution of the underground medium. Multi-channel transient electromagnetic decaying voltage data is input into a pre-trained spatiotemporal adaptive decoupling network. Under the Monte Carlo random dropout mechanism, 100 forward inferences are performed, resulting in 100 sets of three-dimensional predicted resistivity distributions. At the boundary of high-resistivity layers at a depth of 200-300 meters, the transient electromagnetic signal exhibits high nonlinearity when passing through the high-resistivity medium, leading to a large spatial variance value in the 100 prediction outputs. The calculated spatial variation confidence value is only 0.25, lower than the pre-set confidence threshold of 0.5, thus the 200-300 meter depth range is determined to be a prediction blind zone. In the non-uniform adaptive grid construction, the initial grid with a side length of 5 meters is refined into a dense grid with a side length of 2.5 meters in the 200-300 meter depth region using octree partitioning. For shallower regions such as 100 meters, where the confidence level is greater than 0.8, the grid is merged into a sparse grid with a side length of 10 meters.

[0071] The predicted mean resistivity is interpolated and projected onto a non-uniform adaptive grid to form an initial resistivity model. This initial resistivity model predicts a resistivity of 350 Ω·m in the spatial depth region of 200-300 meters, while the actual resistivity of the formation is 500 Ω·m. A spatial weighting factor of 0.25 is correspondingly added to the diagonal weight matrix, reducing the constraint strength of the adaptive regularization matrix in the prediction blind zone of 200-300 meters. Simultaneously, the adjusted damping factor matrix increases the adaptive damping coefficient in the prediction blind zone, limiting drastic fluctuations in model parameter updates. The forward-modeled composite attenuation voltage curve corresponding to the initial resistivity model deviates significantly from the actual observed voltage data. The sensitivity matrix is ​​solved using the adjoint state method under multi-source constraints, and the initial resistivity model and the adaptive regularization matrix are substituted into the Levenberg-McQuart iterative update equation to calculate the model update amount. After 30 iterations, the resistivity value in the region with a space depth of 200 to 300 meters was adjusted from the initial 350 Ω·m to a value close to the actual 500 Ω·m, and the corresponding forward voltage decay curve perfectly matches the actual observation data.

[0072] Reference Figure 2 and Figure 3 , Figure 2This paper compares the decay voltage curves obtained by different methods with the actual observed data. The decay voltage curve obtained by the forward modeling of the data-driven method shows significant discrepancies with the actual observed voltage curve during the observation period of 10 to 40 milliseconds. Particularly in the late decay stage of 30 to 40 milliseconds, the black curve representing the data-driven method deviates significantly from the blue curve representing the actual observed data. This deviation indicates that the three-dimensional resistivity spatial distribution output by the data-driven method does not conform to the physical governing equations of electromagnetic induction, resulting in a lack of theoretical interpretability in the numerical inversion prediction. In contrast, this proposed method introduces a physically constrained objective function for Levenberg-McQuart iteration, ultimately achieving a convergence of the decay voltage curve obtained by the forward modeling with the actual observed data throughout the entire time period.

[0073] Figure 3 This demonstrates a comparison between the conventional Levenberg-McQuart algorithm, based on the assumption of a uniform semi-finite space as the initial model, and the actual values. When using the global average as the initial model of a uniform half-space, the conventional Levenberg-McQuart iterative inversion is prone to getting trapped in local minima due to the strong nonlinearity and multiple solutions inherent in transient electromagnetic inversion. For example... Figure 3 As shown in the four typical resistivity profiles (a, b, c, and d), the blue broken line representing conventional physical inversion exhibits numerous irregular, jagged, and violent oscillations within a depth range of 0 to 1000 meters, completely failing to reconstruct the true stepped resistivity layering structure underground. Especially... Figure 3 In the mid-section (c), conventional physical inversion even produced a resistivity maxima anomaly with a resistivity exceeding 6000 Ω·m near a depth of 600 meters, severely interfering with the boundary localization of deep low-resistivity anomalies. In contrast, this scheme uses the prediction results output by a spatiotemporally adaptive decoupling network as the initial model and constructs a non-uniform adaptive grid using spatial variability confidence, dynamically controlling the damping factor matrix during the iteration process. This control strongly constrains the inversion solution space, effectively eliminating the numerical divergence and multiple-solution oscillations of conventional iterative methods. It not only perfectly reconstructs the true stepped stratum resistivity boundary but also completely eliminates false anomalies at depth, demonstrating the significant technical advantages of this scheme in improving the resolution of shallow and deep boundary identification.

[0074] The present invention also discloses a transient electromagnetic inversion system with multi-scale physical confidence constraints, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the transient electromagnetic inversion method with multi-scale physical confidence constraints as described above.

[0075] The processor can be a central processing unit (CPU). Of course, depending on the actual use, it can also be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), off-the-shelf programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or any conventional processor, etc., and this application does not limit it.

[0076] The memory can be an internal storage unit of a computer device, such as a hard disk or RAM, or an external storage device, such as a plug-in hard disk, smart memory card (SMC), secure digital card (SD), or flash memory card (FC) provided on the computer device. Furthermore, the memory can be a combination of internal storage units and external storage devices of a computer device. The memory is used to store computer programs and other programs and data required by the computer device. The memory can also be used to temporarily store data that has been output or will be output. This application does not limit this.

[0077] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of protection of this application is limited to these examples; within the framework of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of one or more embodiments of this application as described above, which are not provided in detail for the sake of brevity.

[0078] One or more embodiments in this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of this application. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of one or more embodiments in this application should be included within the protection scope of this application.

Claims

1. A transient electromagnetic inversion method with multi-scale physical confidence constraints, characterized in that, The steps include the following: Multi-channel transient electromagnetic data of underground exploration targets excited by secondary induced electromagnetic fields are obtained by electromagnetic sensors. A spatiotemporal adaptive decoupling network containing a time-series encoder and a spatial decoder is constructed. The time-series encoder adopts a cascaded structure of multi-scale temporal convolution and bidirectional long short-term memory network, and the spatial decoder adopts a spatial transposed convolution structure with an attention mechanism. Multi-channel transient electromagnetic data is input into a spatiotemporal adaptive decoupling network, and multiple forward inferences are performed in the spatiotemporal adaptive decoupling network based on the Monte Carlo random drop mechanism to obtain multiple sets of predicted resistivity distributions. The spatial variation confidence level, which represents the uncertainty of the geomorphological prediction of the underground medium where the underground exploration target is located, is obtained by quantifying the spatial distribution variance of the multiple sets of predicted resistivity distributions. Based on the spatial variation confidence level, a non-uniform adaptive grid is constructed for the underground medium where the underground exploration target is located. Multiple sets of predicted resistivity distributions are discretized in the non-uniform adaptive grid to construct an initial resistivity model. The spatial variation confidence is transformed into a spatial weight factor, and the spatial weight factor is dynamically embedded into the damping factor and prior regularization matrix of the LM algorithm in the physical iterative inversion process to construct an objective function containing physical constraints. The adjoint state method under multiple field source constraints is used to solve the sensitivity matrix of the initial resistivity model in a non-uniform adaptive grid, and the objective function is solved through the sensitivity matrix to obtain the model update amount; The initial resistivity model is iteratively updated using model update increments until the preset convergence condition is met, and a three-dimensional resistivity image of the spatial physical boundary of the underground exploration target is output.

2. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 1, characterized in that, Before inputting the multi-channel transient electromagnetic data into the spatiotemporal adaptive decoupling network, the process includes calculating the vertical magnetic field components of a pre-defined subsurface medium model using transient electromagnetic forward modeling equations as simulated observation data, and then using the simulated observation data to construct a training dataset to train the spatiotemporal adaptive decoupling network. The calculation of the vertical magnetic field components using transient electromagnetic forward modeling equations is a key step. The calculation formula is as follows: ， in, Indicates a dimensionless time parameter. Represents the characteristic time constant. Represents the permeability of free space. Indicates the radius of the transmission loop. Indicates the transmitting current. Represents the resistivity of the subsurface medium model. This indicates the preset decay time. This represents the preset error function.

3. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 1, characterized in that, The construction and output process of the time series encoder includes: Multiple preset temporal convolutional kernels with different receptive fields are configured to construct a multi-scale temporal feature extraction layer; Local multi-scale temporal features of multi-channel transient electromagnetic data are extracted using a multi-scale temporal feature extraction layer. A bidirectional long short-term memory network is constructed to capture the temporal dependencies of local multi-scale temporal features in both the forward and reverse temporal directions, thereby obtaining a global temporal feature vector. Local multi-scale temporal features are concatenated and fused with global temporal feature vectors to obtain concatenated temporal features, which are then used as the output of the time series encoder.

4. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 3, characterized in that, The construction and output process of the spatial decoder includes: Construct multiple spatially transposed convolutional layers to progressively upsample the cascaded temporal features to restore them to three-dimensional spatial features; A dual attention unit consisting of channel attention and spatial attention mechanisms is embedded between adjacent spatially transposed convolutional layers; The attention feature weights of the three-dimensional spatial features in the channel dimension and spatial dimension are calculated using dual attention units; Attention feature weights are used to weight and adjust the three-dimensional spatial features to focus on the boundary features of the resistivity spatial distribution, and multiple sets of predicted resistivity distributions are output.

5. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 1, characterized in that, The method of performing multiple forward inferences in a spatiotemporally adaptive decoupled network based on the Monte Carlo random dropout mechanism to obtain multiple sets of predicted resistivity distributions, and obtaining the spatial variation confidence level representing the uncertainty of the geomorphological prediction of the subsurface medium where the subsurface exploration target is located by statistically quantifying the spatial distribution variance of multiple sets of predicted resistivity distributions, includes the following steps: Configure Monte Carlo random dropout layers in the hidden layers of the spatiotemporal adaptive decoupling network; During the network prediction phase, the Monte Carlo random drop-out layer is kept active. Multi-channel transient electromagnetic data is repeatedly input into the spatiotemporal adaptive decoupling network to perform forward prediction inference a preset number of times, resulting in multiple sets of predicted resistivity distributions. Calculate the average resistivity and spatial distribution variance of multiple predicted resistivity distributions at each spatial node; The reciprocal of the spatial distribution variance is determined as the spatial variation confidence level, representing the uncertainty of predicting the geomorphology of the underground medium where the underground exploration target is located.

6. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 1, characterized in that, The process of constructing a non-uniform adaptive grid for the subsurface medium where the subsurface exploration target is located based on spatial variability confidence, and discretizing multiple sets of predicted resistivity distributions within the non-uniform adaptive grid to construct an initial resistivity model includes the following steps: Obtain a preset background spatial grid containing the initial subdivision accuracy, and match the corresponding spatial variation confidence at the nodes of the background spatial grid; In the background spatial grid, identify prediction blind zones where the confidence level of spatial variation is less than a preset confidence threshold, and increase the initial subdivision accuracy within the prediction blind zones; In the background spatial grid, identify confidence regions with spatial variation confidence not less than the confidence threshold, and reduce the initial subdivision accuracy within the confidence regions to construct a non-uniform adaptive grid for the subsurface medium where the subsurface exploration target is located. Multiple sets of predicted resistivity distributions are interpolated and projected onto a non-uniform adaptive grid for discretization to obtain the initial resistivity model.

7. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 6, characterized in that, The step of converting spatial variation confidence into spatial weighting factors and dynamically embedding these spatial weighting factors into the damping factor and prior regularization matrix of the LM algorithm during the physical iterative inversion process to construct an objective function containing physical constraints includes the following steps: The spatial variation confidence level at each grid cell in the initial resistivity model is determined as the spatial weighting factor. Construct a diagonal matrix to constrain the spatial inversion weights using spatial weighting factors; By using a diagonal matrix to perform a product weighting adjustment on a pre-defined prior regularization matrix, the regularization constraint strength of a non-uniform distribution can be obtained. The preset damping factor matrix of the LM algorithm during the iteration process is adjusted by using a diagonal matrix to increase the adaptive damping coefficient at the prediction blind zone, so as to obtain the adjusted damping factor matrix. The regularization constraint strength and the adjusted damping factor matrix are used as physical confidence constraints and embedded in the physical iteration of the LM algorithm to construct the objective function.

8. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 7, characterized in that, The objective function formula is as follows: ， in, This represents the value of the objective function. This represents the current model parameters after discretization on a non-uniform adaptive grid. This represents multi-channel transient electromagnetic data. Indicates the current model parameters Forward modeling composite data on a non-uniform adaptive grid Represents the adaptive regularization matrix. This represents the preset regularization parameters.

9. The transient electromagnetic inversion method with multi-scale physical confidence constraints according to claim 8, characterized in that, The specific model parameter update equation corresponding to the objective function is as follows: ， in, Indicates the first The sensitivity matrix corresponding to the next iteration Indicates the first The amount of model parameter updates in each iteration. Indicates the first Model parameters for the next iteration Indicates the first The forward-modeled composite data of the model parameters in the next iteration on a non-uniform adaptive grid. This represents the preset damping control coefficient. This represents the adjusted damping factor matrix.

10. A transient electromagnetic inversion system with multi-scale physical confidence constraints, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the transient electromagnetic inversion method with multi-scale physical confidence constraints as described in any one of claims 1 to 9.

Images