A two-dimensional magnetotelluric inversion method

Abstract

This application relates to the technical field of exploration and data processing, and in particular to a magnetotelluric two-dimensional inversion method. The method includes: acquiring and normalizing magnetotelluric measured impedance data to construct network input data; constructing a GU-CBAM-Unet inversion network model, which is based on the Unet architecture, replacing batch normalization layers with group normalization GN layers to adapt to single-sample training scenarios, and introducing CBAM attention modules at downsampling paths and skip connections to enhance anomaly feature extraction capabilities; constructing a total loss function containing fitting error terms and model roughness constraint terms based on a differentiable magnetotelluric two-dimensional forward modeling operator; inputting single-sample measured data into the network and driving unsupervised iterative training with the total loss function; and outputting a high-resolution two-dimensional resistivity model after training convergence. This application does not require massive amounts of labeled data; it can achieve stable and high-resolution subsurface two-dimensional resistivity model inversion using only measured data from a single exploration.

Description

Technical Field

[0001] This application relates to the technical field of exploration and data processing, and in particular to a two-dimensional magnetotelluric inversion method. Background Technology

[0002] Magnetotelluric (MT) sounding is an important geophysical exploration method for probing the deep electrical structure of the Earth and searching for deep mineral resources. In recent years, deep learning technology, due to its powerful nonlinear mapping capabilities, has been introduced into two-dimensional magnetotelluric inversion.

[0003] Currently, existing deep learning-based magnetotelluric inversion methods (such as those based on conventional Unet networks) typically employ a supervised learning paradigm. This method requires pre-constructing a massive paired dataset containing corresponding "theoretical electromagnetic responses" and "actual subsurface resistivity distributions" as "labels" to pre-train the neural network. After the network converges, measured magnetotelluric data is then input into the network to predict the resistivity model.

[0004] However, in actual geophysical exploration, it is impossible to obtain massive amounts of real subsurface geological samples (labels) to train the network. When attempting to directly use unlabeled measured data from a single exploration for unsupervised iterative optimization, conventional networks will completely fail due to the extremely small sample size. Specifically, conventional deep learning networks (such as UNet) generally contain batch normalization layers (BatchNorm, BN). When only a single sample or a very small batch of data (e.g., measured data from a single survey line) is input into the network, the BN layer's estimation of statistical quantities such as the mean and variance of the data becomes extremely unstable and inaccurate, leading to chaotic gradients and failure to converge during the iteration process. At the same time, due to the lack of reasonable physical constraints, conventional networks cannot spontaneously infer reasonable subsurface structures from single sample data, ultimately resulting in the failure of high-resolution inversion.

[0005] Therefore, existing neural networks rely on massive amounts of labeled data and suffer from single-sample training failure due to their internal batch normalization mechanism. As a result, they cannot achieve stable and high-resolution underground two-dimensional resistivity model inversion using unlabeled measured data (single sample) obtained from a single exploration. Summary of the Invention

[0006] To address the aforementioned problems, this application provides a two-dimensional magnetotelluric inversion method, employing the following technical solution: Includes the following steps: Step S1: Obtain the measured magnetotelluric impedance data of the exploration area to be inverted, and normalize the measured impedance data to construct the network input data; at the same time, obtain the coordinates and frequency points of the resistivity model subdivision nodes of the exploration area. Step S2: Construct a GU-CBAM-Unet inversion network model; the GU-CBAM-Unet inversion network model adopts the Unet basic network structure, replaces the batch normalization layer in the basic network structure with the group normalization GN layer, and introduces the convolutional attention CBAM module at the downsampling path and skip connection; Step S3: Based on the differentiable magnetotelluric two-dimensional forward modeling operator, construct a total loss function that includes fitting error terms and model roughness constraint terms; Step S4: Input the network input data, resistivity model partitioning node coordinates and frequency point data as single samples into the GU-CBAM-Unet inversion network model, and drive the network to perform iterative training with the total loss function to optimize the weight parameters of the network model; Step S5: Stop training when the total loss function converges, and output the feature data of the finally converged GU-CBAM-Unet inversion network model as a high-resolution two-dimensional resistivity model of the target geological exploration area.

[0007] By adopting the above technical solution, this method uses the GU-CBAM-Unet inversion network model. By replacing the batch normalization layer with the group normalization (GN) layer, the dependence on batch dimension statistics is eliminated. This allows the network to maintain the stability of feature distribution and effective gradient propagation even in unsupervised training scenarios with only single-sample measured data. This fundamentally solves the training collapse problem caused by the inaccuracy of BN layer statistics in conventional networks under extremely small sample conditions. At the same time, the CBAM attention module is introduced to enhance the network's ability to focus on key geophysical features. Iterative optimization is driven by a physical constraint loss function based on differentiable forward modeling operators. Under conditions where no label data is required, stable convergence and output of a high-resolution two-dimensional resistivity model can be achieved using only measured data from a single exploration. Optionally, in step S1, the magnetotelluric measured impedance data includes apparent resistivity data and impedance phase data in TM polarization mode and / or TE polarization mode.

[0008] By adopting the above technical solution and using apparent resistivity and impedance phase data in TM and / or TE polarization modes as input, complementary information on underground electrical structures can be obtained from different polarization angles. The TM mode is sensitive to vertical electrical interfaces while the TE mode is sensitive to horizontal electrical interfaces. The combined use of the two modes can more comprehensively constrain the inversion model and effectively improve the resolution of the boundaries and depths of underground anomalies in the inversion results.

[0009] Optionally, in step S1, the measured impedance data is normalized, specifically by using the formula: Normalization is performed; for the impedance phase data, the formula is: Normalization is performed; among which, The normalized apparent resistivity To measure the apparent resistivity, and These are the maximum and minimum values ​​of the measured logarithmic resistivity, respectively. The normalized impedance phase, The measured impedance phase, and These represent the maximum and minimum values ​​of the measured impedance phase, respectively.

[0010] By adopting the above technical solutions, logarithmic domain normalization is used for apparent resistivity because the numerical variation range of underground resistivity usually spans several orders of magnitude. Directly using linear values ​​would lead to large numerical data dominating the gradient during network training. Taking the logarithm and then performing max-min normalization can compress apparent resistivity data of different orders of magnitude into a uniform numerical range. Max-min normalization is directly performed for impedance phase because phase data usually varies within a fixed range and there is no problem of crossing orders of magnitude. Both normalization methods map the data to the [0,1] interval, which can accelerate network convergence and improve training stability.

[0011] Optionally, when using TM polarization mode combined with TE polarization mode inversion, the steps for constructing the network input data include: sequentially merging and splicing the normalized TM apparent resistivity, TM impedance phase, TE apparent resistivity and TE impedance phase to form a four-channel multidimensional tensor with a shape of 1×4×number of frequency points×number of observation points, which serves as the single-sample training input data for the network.

[0012] By adopting the above technical solution, the four normalized data are constructed as a four-channel tensor input, enabling the network to simultaneously receive the apparent resistivity and phase information of two polarization modes in a single forward propagation. This fully utilizes multi-channel convolution operations to automatically learn the correlation features and complementary information between data of different polarization modes, avoiding the information fragmentation problem when processing multiple types of data separately, thereby improving the network's comprehensive discrimination ability for complex underground electrical structures.

[0013] Optionally, in step S2, the group normalization GN layer does not rely on the batch dimension of the input data for statistics. Instead, it groups the channel features and calculates the mean and variance within each group for normalization to ensure the stability of the feature distribution of single-sample measured data during training iterations.

[0014] By adopting the above technical solution, the group normalization (GN) layer groups the channel features and calculates the normalized statistics independently within each group. Its normalization process does not depend on the number of samples in the batch dimension. Therefore, even in a single-sample training scenario with a batch size of 1, the mean and variance within each channel group can still be reliably estimated, thereby ensuring the stability of the feature distribution in the intermediate layers of the network. This avoids the problem of gradient chaos and training collapse caused by drastic fluctuations in statistics in conventional BN layers under single-sample conditions, making unsupervised single-sample iterative optimization possible.

[0015] Optionally, in step S2, the convolutional attention CBAM module is a dual-path structure, consisting of a cascaded channel attention submodule and a spatial attention submodule; this module is used to suppress the geophysical background interference features extracted by the network and focus on the spatial geometric location features of deep electrical anomalies.

[0016] By adopting the above technical solution, the CBAM module adaptively weights the features of different channels through the channel attention submodule, suppressing channels containing redundant background information and enhancing channels containing anomaly information. Then, through the spatial attention submodule, it adaptively emphasizes the feature response of the region where the anomaly is located in the spatial dimension. The two submodules work together in series, which can effectively guide the network to focus on the spatial geometric location and boundary contour features of deep electrical anomalies, reduce the interference of the uniform background medium on the inversion accuracy, and thus significantly improve the resolution of the inversion results and the accuracy of anomaly localization.

[0017] Optionally, the internal processing logic of the Convolutional Attention (CBAM) module is executed according to the following formula: F′=M c (F)⨂F(1) F′′=M s (F′)⨂F′(2) M c (F)= σ ( MLP ( AvgPool (F))+ MLP ( MaxPool (F)))(3) M s (F)= σ ( f 7×7 ([ AvgPool (F); MaxPool (F)]))(4) Wherein, the input feature map F∈R C×H×W F′ is the feature map after channel attention weighting, and M c ∈R C×1×1 For channel attention maps, M s ∈R1×H×W F′′ is the spatial attention map, and F′′ is the final output feature map. σ It is the Sigmoid activation function. f 7 ×7 For 7×7 size convolution operations, MLP This represents a multilayer perceptron. AvgPool and MaxPool ∠ and ∠ represent average pooling and max pooling operations, respectively, and ⨂ represents element-wise multiplication.

[0018] By adopting the above technical solutions, channel attention captures the global average response and peak response information in the channel dimension by simultaneously using average pooling and max pooling, and then fuses them to generate channel weights after MLP mapping, ensuring a comprehensive evaluation of the importance of different feature channels. Spatial attention compresses the feature map after channel attention weighting by average pooling and max pooling in the spatial dimension, and then fuses spatial context information through a 7×7 large-size convolution to generate a spatial weight map. The use of large-size convolution kernels expands the spatial receptive field, enabling the attention mechanism to capture a wider range of spatially related features. The two-stage element-wise multiplication operation realizes fine-grained adaptive calibration of the feature map.

[0019] Optionally, in step S3, the total loss function is calculated using the following formula: ;in, Total loss; This represents the fitting error term between the apparent resistivity forward response of the resistivity model output by the network and the measured normalized apparent resistivity. The term representing the fitting error between the impedance phase forward response of the resistivity model output by the network and the measured normalized impedance phase; m∈{TM,TE} represents the polarization mode; Θ represents the roughness constraint term of the model; Θ represents the two-dimensional resistivity model output by the network; ∇x and ∇z represent the gradients in the transverse and depth directions of the survey line, respectively; This indicates the calculation of the L2 norm; λ is the model smoothness weighting coefficient determined based on the environmental noise level.

[0020] By adopting the above technical solutions, the data fitting error term in the total loss function is converted into a theoretical impedance response by a differentiable forward modeling operator and compared with the measured data. This realizes the direct driving force of physical laws on network training, enabling the network to spontaneously search for a reasonable resistivity model by minimizing the deviation between the forward response and the measured data without label data. The model roughness constraint term suppresses unreasonable violent oscillations and spurious anomalies in the inversion solution by applying L2 norm penalties to the output model in the lateral and depth directions, guiding the model to tend towards smoothness and continuity, which conforms to the natural laws of the electrical distribution of geological bodies. The introduction of the λ coefficient allows for flexible adjustment of the smoothness degree according to the noise level of the measured data, achieving a balance between data fitting and model rationality.

[0021] Optionally, the fitting error term is calculated using the mean absolute error (MAE), and the specific formula is as follows: ; ;in, The differentiable magnetotelluric two-dimensional forward modeling operator is used to convert the two-dimensional resistivity model Θ output by the current iterative network into the corresponding theoretical impedance response in real time. and Indicates the measured magnetotelluric apparent resistivity and impedance phase; This indicates a normalization operation. This indicates the calculation of the mean absolute error.

[0022] By adopting the above technical solution, the theoretical response is calculated in real time using a differentiable forward operator, ensuring that the gradient can be backpropagated from the loss function to the network weights via forward calculation, thus achieving end-to-end optimization. The forward theoretical response is normalized in the same way as the measured data before calculating MAE, ensuring that the two are compared on the same numerical scale and avoiding error bias caused by inconsistency in dimensions. MAE is selected as the error metric, which is more robust to outliers than mean squared error, and can reduce the interference of a few abnormal noise points in the measured data on the overall training process, thereby improving the robustness of inversion.

[0023] In summary, this application includes at least one of the following beneficial technical effects: 1. By replacing the batch normalization layer in the Unet network with the group normalization GN layer, the network's dependence on batch dimension statistics is fundamentally eliminated. This allows the network to maintain stable feature distribution and effective gradient propagation even under training conditions with only single-sample unlabeled measured data. Combined with the physical constraint loss function constructed by the differentiable forward model to drive unsupervised iterative optimization, the goal of stably reproducing a high-resolution two-dimensional resistivity model using only single exploration measured data without the need for massive amounts of labeled data is achieved. 2. By introducing the CBAM attention module at the downsampling path and skip connection, the network can adaptively filter important features in the channel dimension and focus on the anomalous body region in the spatial dimension, effectively suppressing the interference features of geophysical background and significantly improving the spatial positioning accuracy and boundary resolution of deep electrical anomalies in the inversion results. Attached Figure Description

[0024] Figure 1 This is a schematic diagram of the overall process of the magnetotelluric two-dimensional inversion method according to an embodiment of this application; Figure 2 This is a schematic diagram of the architecture of the GU-CBAM-Unet inversion network model according to an embodiment of this application; Figure 3 This is a schematic diagram of the internal structure of the CBAM module according to an embodiment of this application; Figure 4 This is a schematic diagram of unsupervised iterative training driven by the total loss function in an embodiment of this application. Detailed Implementation

[0025] The following is in conjunction with the appendix Figure 1 - Appendix Figure 4 This application will be described in further detail below.

[0026] A two-dimensional magnetotelluric inversion method, referring to Figure 1 The method includes the following steps S1 to S5: Step S1: Obtain the measured data and perform preprocessing.

[0027] First, magnetotelluric impedance data of the exploration area to be inverted is acquired. In this embodiment, the magnetotelluric impedance data includes apparent resistivity data and impedance phase data under TM polarization mode (transverse magnetic mode) and TE polarization mode (transverse electric mode). Specifically, a magnetotelluric sounding station array is deployed in the field, and natural electromagnetic field signals are recorded using magnetotelluric sounding instruments. After time series processing and impedance tensor estimation, the apparent resistivity of each observation station at different frequencies in TM mode is obtained. TM mode impedance phase TE mode apparent resistivity and TE mode impedance phase .

[0028] Assuming that Ns observation stations are set up along the survey line in the exploration area, and each station collects data at Nf frequency points, then the apparent resistivity and impedance phase data under each polarization mode constitute an Nf×Ns two-dimensional matrix.

[0029] Simultaneously, the resistivity model partitioning node coordinates and frequency data of the exploration area are acquired. The resistivity model partitioning node coordinates refer to the node coordinates of the two-dimensional discrete grid used for finite difference forward modeling, including the node coordinate sequences along the survey line direction (x-direction) and the depth direction (z-direction). This grid defines the spatial discretization scheme of the two-dimensional resistivity model to be inverted. The frequency data are the frequency sequences corresponding to the measured data. This information will be used in subsequent steps for calculating the differentiable forward modeling operator.

[0030] Then, the measured impedance data is normalized. The purpose of normalization is to map the apparent resistivity and impedance phase data, which have large differences in numerical range, to the [0,1] interval, so as to improve the stability and convergence efficiency of network training.

[0031] Specifically, for apparent resistivity data, since subsurface resistivity typically spans several orders of magnitude (e.g., from 1 Ω·m to 10000 Ω·m), directly using the raw values ​​would lead to numerical scaling imbalance. Therefore, the apparent resistivity is first taken as a base-10 logarithm, and then normalized using a min-max normalization. The normalization formula is:

[0032] in, The normalized apparent resistivity To measure the apparent resistivity, and These are the maximum and minimum values ​​of the logarithmic apparent resistivity among all measured values ​​under this polarization mode.

[0033] For impedance phase data, since the phase typically varies within the range of 0° to 90° (for TM mode) or within a certain finite range (for TE mode), there is no issue of crossing orders of magnitude; therefore, maximum-minimum normalization is directly applied. The normalization formula is:

[0034] in, The normalized impedance phase, The measured impedance phase, and These are the maximum and minimum values ​​among all measured impedance phases in this polarization mode, respectively.

[0035] After normalization, the normalized data is used to construct the network input data. In this embodiment, the TM polarization mode is combined with the TE polarization mode for inversion. The steps for constructing the network input data are as follows: the normalized TM apparent resistivity matrix (Nf×Ns), TM impedance phase matrix (Nf×Ns), TE apparent resistivity matrix (Nf×Ns), and TE impedance phase matrix (Nf×Ns) are sequentially merged and concatenated along the channel dimension to form a four-channel multidimensional tensor with a shape of 1×4×Nf×Ns. The batch dimension is 1 (i.e., a single sample), the channel dimension is 4 (corresponding to TM apparent resistivity, TM impedance phase, TE apparent resistivity, and TE impedance phase, respectively), the height dimension is the number of frequency points Nf, and the width dimension is the number of observation points Ns. This four-channel tensor serves as the single-sample training input data for the GU-CBAM-Unet inversion network model.

[0036] It should be noted that when using only a single polarization mode (TM mode or TE mode) for inversion, the number of input tensor channels is reduced to 2 (apparent resistivity channel and impedance phase channel), with a shape of 1×2×Nf×Ns.

[0037] Step S2: Construct the GU-CBAM-Unet inversion network model.

[0038] Reference Figure 2 The GU-CBAM-Unet inversion network model uses the Unet basic network structure as its skeleton and makes two key improvements on this basis: First, all batch normalization (BN) layers in the basic network structure are replaced with group normalization (GN) layers; Second, convolutional block attention (CBAM) modules are introduced at the downsampling path and skip connection.

[0039] Regarding the basic network structure of Unet: The Unet network is a classic fully convolutional network with an encoder-decoder structure. It consists of symmetrical downsampling paths (encoder) and upsampling paths (decoder), and uses skip connections to concatenate and fuse the feature maps of the encoder layers with the feature maps of the corresponding decoder layers.

[0040] The encoder consists of multiple downsampling stages. Each downsampling stage contains two convolutional blocks (each convolutional block consists of a convolutional layer, a normalization layer, and an activation function layer in sequence) and a downsampling operation (such as max pooling or a convolution with a stride of 2). In this embodiment, the encoder contains four downsampling stages. The convolutional kernel size is uniformly set to 3×3, and the padding is set to 1 to maintain the feature map space size. The number of output channels in each stage is 64, 128, 256, and 512, respectively. At the bottom layer (bottleneck layer) of the encoder, a convolutional block with 1024 output channels is set.

[0041] The decoder also consists of multiple upsampling stages. Each upsampling stage includes an upsampling operation (such as transposed convolution or bilinear interpolation upsampling followed by convolution), which doubles the spatial size of the feature map. This feature map is then concatenated with the feature map from the corresponding encoder layer via skip connections in the channel dimension, and finally fused and processed by two convolutional blocks. The number of output channels for each upsampling stage is 512, 256, 128, and 64, respectively.

[0042] At the end of the decoder, a 1×1 convolutional layer is set as the output layer, which maps the number of channels to 1, that is, outputs a single-channel two-dimensional feature map. The spatial size of the feature map is consistent with the preset resistivity model grid size, and each pixel value represents the logarithmic resistivity value of the corresponding grid cell.

[0043] A sigmoid activation function is applied after the output layer to restrict the output value to the (0,1) interval, and then it is mapped to a preset logarithmic range of resistivity through a linear scaling transformation. Where ρmin and ρmax are the lower and upper limits of the preset resistivity search range, respectively. The specific scaling formula is:

[0044] Where x is the original output of the 1×1 convolutional layer. This represents the logarithmic value of the final output resistivity. This design ensures that the resistivity value output by the network is always within a physically reasonable range.

[0045] It should be noted that since the spatial size of the network input data (Nf×Ns, i.e., number of frequency points × number of observation points) is usually inconsistent with the preset resistivity model grid size (Nz×Nx, i.e., number of grids in the depth direction × number of grids in the survey line direction), appropriate size adjustment operations (such as interpolation scaling or padding clipping) need to be performed at the input end of the encoder or the output end of the decoder to ensure smooth data flow in the network and the matching of the output feature map size with the target model grid size.

[0046] Regarding the replacement of the group-normalized GN layer: In a standard Unet network, each convolutional block typically employs a batch normalization (BN) layer. During training, the BN layer calculates the mean and variance of the current mini-batch along the batch dimension to normalize the feature maps. When the batch size is large (e.g., 32 or 64), the BN layer provides stable and reliable statistical estimates, effectively accelerating training and improving generalization ability.

[0047] However, in the application scenario of this application, the input data is single-sample measured data obtained from a single exploration, and the batch size is fixed at 1. Under this condition, the BN layer has only one sample in the batch dimension, and the estimation of its mean and variance degenerates into the statistics of that single sample, completely losing the significance of the statistical average between batches. This leads to the extremely unstable feature distribution after normalization, which in turn causes drastic fluctuations in gradients or even complete training failure during the network training process.

[0048] To address the aforementioned issues, this application replaces the BN layers within all convolutional blocks in the Unet basic network with group normalized (GN) layers. These GN layers do not rely on the batch dimension of the input data for statistical analysis; instead, they group channel features and calculate the mean and variance within each group for normalization.

[0049] Specifically, for a feature map of shape [B,C,H,W] (where B is the batch size, C is the number of channels, and H and W are the spatial height and width), the GN layer uniformly divides the C channels into G groups, each containing C / G channels. Then, within each group, the mean μg and variance are independently calculated along the channel and spatial dimensions (i.e., over C / G×H×W elements). And normalize it:

[0050] Where xi is the eigenvalue belonging to the g-th group, and ϵ is a very small positive number (e.g., ...). This is used to prevent division by zero. After normalization, a learnable scaling parameter γ and an offset parameter β are also added for affine transformation.

[0051] In this embodiment, the number of groups G is set to 32. Since the statistical calculation of the GN layer is performed entirely on the intra-group channels and spatial dimensions of a single sample, and is independent of the batch size B, when B=1, the GN layer can still obtain stable and reliable statistical estimates based on a sufficient number of elements (C / G×H×W), thereby ensuring the stability of the feature distribution of single-sample measured data during training iterations and enabling the network to converge stably.

[0052] Regarding the introduction of the Convolutional Attention (CBAM) module: Reference Figure 3The Convolutional Attention (CBAM) module has a dual-path structure, consisting of a cascaded Channel Attention Module and a Spatial Attention Module.

[0053] In this application, a CBAM module is introduced into the downsampling path and skip connections of the GU-CBAM-Unet inversion network model. Specifically, a CBAM module is inserted after two convolutional blocks in each downsampling stage of the encoder and before the downsampling operation; simultaneously, a CBAM module is also applied to the feature map before the skip connections pass the encoder feature map to the decoder for concatenation. This module is used to suppress the geophysical background disturbance features extracted by the network and focus on the spatial geometric location features of deep electrical anomalies.

[0054] The internal processing logic of the CBAM module is executed according to the following steps: Step 1: Channel attention calculation.

[0055] Given an input feature map F∈R C×H×W (Where C is the number of channels, and H and W are the spatial height and width of the feature map), first, global average pooling (AvgPool) and global max pooling (MaxPool) are performed on F respectively to compress the spatial information of each channel into a scalar value, resulting in two channel description vectors. and .

[0056] These two description vectors are then input into a multilayer perceptron (MLP) with shared weights. The MLP consists of two fully connected layers: the first layer compresses the number of channels from C to C / r (where r is the compression ratio, which is set to 16 in this embodiment) and applies the ReLU activation function; the second layer restores the number of channels from C / r to C.

[0057] Finally, the outputs of the MLP on the two description vectors are summed element-wise, and then passed through a sigmoid activation function to generate a channel attention map M. c ∈R C×1×1 Its mathematical expression is: M c (F)= σ ( MLP ( AvgPool (F))+ MLP ( MaxPool (F))) Each element in the channel attention map Mc(F) has a value between 0 and 1, representing the importance of the corresponding channel feature.

[0058] The channel attention map is then multiplied element-wise with the input feature map (the channel attention map is automatically broadcast along the spatial dimensions) to obtain the channel attention-weighted feature map: F′=M c (F)⨂F Step 2: Spatial attention calculation.

[0059] The feature map F′∈R after channel attention weighting C×H×W Average pooling and max pooling are performed along the channel dimension to obtain two single-channel spatial description graphs. and .

[0060] The two spatial description maps are concatenated along the channel dimension to form a 2×H×W feature map, and then a convolution operation with a kernel size of 7×7 and padding of 3 is applied to it. f 7×7 (The number of output channels is 1), and then a spatial attention map M is generated by passing it through the Sigmoid activation function. s ∈R 1×H×W Its mathematical expression is: M s (F′)= σ ( f 7×7 ([ AvgPool (F′); MaxPool (F′)])) The square brackets [⋅;⋅] indicate concatenation along the channel dimension. Each element in the spatial attention map has a value between 0 and 1, representing the importance of the corresponding spatial location feature. A large 7×7 convolutional kernel is chosen to obtain a larger receptive field in the spatial dimension, thereby better capturing spatial contextual information.

[0061] Finally, the spatial attention map and the channel attention-weighted feature map are multiplied element-wise (the spatial attention map is automatically broadcast along the channel dimension) to obtain the final output feature map of the CBAM module: F′′=M s (F′)⨂F′ The output feature map Compared to the original input feature map F, the feature channels containing important geophysical information are emphasized in the channel dimension, and the feature responses of the region where the anomaly is located are emphasized in the spatial dimension, thereby enhancing the network's ability to locate and identify deep electrical anomalies.

[0062] Step S3: Construct the total loss function.

[0063] Reference Figure 4This step is based on a differentiable magnetotelluric two-dimensional forward modeling operator to construct a total loss function that includes fitting error terms and model roughness constraint terms, so as to drive the network to perform unsupervised iterative optimization under unlabeled conditions.

[0064] Regarding the differentiable two-dimensional forward magnetotelluric operator: The differentiable magnetotelluric two-dimensional forward modeling operator This is the core physical module for achieving unsupervised inversion in this method. The function of this operator is: given a two-dimensional resistivity model Θ (output by the network), the coordinates of the partition nodes of the resistivity model, and frequency data, it calculates the theoretical electromagnetic response of the resistivity model at each observation station and frequency point by solving the partial differential equation of the two-dimensional forward magnetotelluric problem (Helmholtz equation discretized by the finite difference method), including the theoretical apparent resistivity and theoretical impedance phase.

[0065] The key to this forward operator lies in its "differentiability." In its implementation, the forward operator is constructed as part of an automatically differentiable computation graph within a deep learning framework (such as PyTorch). Specifically, all the mathematical operations involved in the forward calculation—assembling finite difference equations, solving linear equations, and calculating impedance, apparent resistivity, and phase from electromagnetic field components—are implemented using tensor operations and automatic differentiation mechanisms provided by the deep learning framework. Thus, when backpropagating the network with the total loss function as the objective function, the gradient can be transferred from the loss function via the computation of the forward operator. Figure 1 The gradient is propagated back to the network's weight parameters, thus achieving end-to-end gradient optimization.

[0066] In each training iteration, the network outputs a two-dimensional resistivity model Θ, and a forward modeling operator. The model is converted into theoretical apparent resistivity and theoretical impedance phase under the corresponding polarization mode in real time. Then, the theoretical response data obtained by forward modeling is normalized in the same way as the measured data, and the error is calculated by comparing it with the normalized measured data.

[0067] Regarding the fitting error term: The fitting error term measures the deviation between the forward-modeled theoretical response of the resistivity model currently output by the network and the measured normalized data. In this embodiment, the fitting error term is calculated using the mean absolute error (MAE), and the specific formula is as follows:

[0068]

[0069] in, The differentiable magnetotelluric two-dimensional forward modeling operator is used to convert the two-dimensional resistivity model Θ output by the current iterative network into the corresponding theoretical impedance response in real time. and Indicates the measured magnetotelluric apparent resistivity and impedance phase; This indicates a normalization operation. This indicates the calculation of the mean absolute error.

[0070] The reason for choosing MAE instead of MSE is that MAE is more robust to outliers in the data. Inevitably, some outliers severely affected by noise exist in measured magnetotelluric data. MSE amplifies the error of these outliers by square, causing the training process to be dominated by a few outliers. MAE, on the other hand, performs an equal-weighted linear summation of the errors of all data points, effectively reducing the interference of outliers on the overall optimization direction.

[0071] Regarding model roughness constraints: Magnetotelluric inversion is essentially an ill-posed problem, meaning that there are infinitely many resistivity models that can fit the same set of measured data within a given error range. To select the most reasonable model from these equivalent solutions, prior constraints need to be introduced. In this embodiment, model roughness constraints (also known as smoothness regularization) are used. The physical assumption is that the actual subsurface resistivity distribution typically exhibits a certain spatial continuity and should not show unreasonable, drastic jumps.

[0072] The formula for calculating the model roughness constraint term is:

[0073] Where Θ represents the two-dimensional resistivity model output by the network; ∇x and ∇z represent the gradients in the transverse and depth directions of the survey line, respectively; This indicates the calculation of the L2 norm; λ is the model smoothness weighting coefficient determined based on the environmental noise level. In actual calculations, It is approximated by the difference between adjacent grid cells along the x-direction of Θ. It is approximated by the difference between adjacent grid cells along the z-direction of Θ.

[0074] λ is the model smoothness weight coefficient, also known as the regularization parameter, and its value is determined based on the environmental noise level of the exploration area. When the noise level of the measured data is high, the value of λ should be appropriately increased to enhance the smoothing constraint and prevent the network from introducing spurious anomalies into the model due to overfitting to noisy data. When the quality of the measured data is good and the noise level is low, the value of λ can be appropriately decreased to allow the model to retain more detailed features. In this embodiment, the typical range of λ values ​​is... The specific value is determined through experiments based on the signal-to-noise ratio of actual exploration data.

[0075] Regarding the complete form of the total loss function: Combining the data fitting error term and the model roughness constraint term, we obtain the complete expression for the total loss function:

[0076] in, Total loss; This represents the fitting error term between the apparent resistivity forward response of the resistivity model output by the network and the measured normalized apparent resistivity. This represents the fitting error term between the impedance phase forward modeling response of the resistivity model output by the network and the measured normalized impedance phase; m∈{TM,TE} represents the polarization mode. When only a single polarization mode is used, the total loss is reduced to three parts: the two fitting error terms of that mode plus the model roughness constraint term.

[0077] The design concept of the total loss function is: the data fitting error term ensures that the resistivity model output by the network can physically explain the measured electromagnetic response data, and the model roughness constraint term ensures that the output model has reasonable smoothness in space. The combined driving force of the two enables the network to spontaneously search for a two-dimensional resistivity model that can fit the measured data and has geological rationality without the need for labels.

[0078] Step S4: Iterative training and optimization.

[0079] The network input data (a four-channel tensor of 1×4×Nf×Ns) constructed in step S1, the node coordinates and frequency data of the resistivity model partitioning are used as single samples to input the GU-CBAM-Unet inversion network model constructed in step S2. The total loss function constructed in step S3 drives the network to perform iterative training in order to optimize the weight parameters of the network model.

[0080] In this embodiment, the Adaptive Moment Estimation (Adam) optimizer is used to iteratively optimize the weight parameters of the GU-CBAM-Unet inversion network model. The Adam optimizer adaptively adjusts the learning rate of each parameter by combining the first moment (mean) and second moment (uncentered variance) of the gradient. Its update rule is as follows:

[0081]

[0082]

[0083]

[0084] Where gt is the gradient at step t, and mt and vt are the estimates of the first and second moments, respectively. and For the bias-corrected moment estimate, β1 and β2 are the exponential decay rates of the moment estimate (usually set to 0.9 and 0.999), η is the learning rate, and ϵ is a very small positive number (e.g., ...). Prevent division by zero.

[0085] During iterative training, a cosine annealing strategy is used to dynamically adjust the learning rate decay. The cosine annealing strategy causes the learning rate to smoothly decay from the initial learning rate ηmax to the minimum learning rate ηmin throughout the training process, following the shape of a cosine function. The formula is as follows:

[0086] Where ηt is the learning rate at iteration t, and T is the total number of iterations (or the length of one cosine annealing cycle). In this embodiment, the initial learning rate ηmax is set to... The minimum learning rate ηmin is set to The total number of iterations T is set according to the actual convergence situation, with a typical value of 2000 to 10000 steps.

[0087] The specific process for each iteration is as follows: (1) Input the network input data into the GU-CBAM-Unet inversion network model, perform forward propagation, and the network outputs a two-dimensional resistivity model Θ; (2) Input the resistivity model Θ along with the coordinates of the subdivision nodes and the frequency data into a differentiable two-dimensional magnetotelluric forward modeling operator. The theoretical apparent resistivity and theoretical impedance phase corresponding to the model are calculated. (3) Normalize the theoretical response obtained by forward modeling, and then calculate the fitting error term between the theoretical response and the measured normalized data; (4) Calculate the roughness constraint term of the resistivity model Θ; (5) Add the fitting error term and the roughness constraint term to obtain the total loss. ; (6) Backpropagate the total loss and calculate the gradient of the total loss with respect to all learnable parameters of the network. During the backpropagation process, the gradient passes through the computation graph of the forward operator and the computation graph of the network in turn to realize end-to-end gradient calculation. (7) The Adam optimizer updates the network's weight parameters based on the calculated gradients; (8) The cosine annealing strategy is used to update the learning rate for the next iteration; (9) Repeat steps (1) to (8) above until the stopping condition is met.

[0088] It is important to emphasize that throughout the entire iterative training process, the same single-sample measured data is used as both the input and the fitting target, without involving any external labeled data. The network's weight parameters are initialized randomly and automatically optimized entirely by the physical constraint loss function; this is the core mechanism of our unsupervised inversion method.

[0089] Step S5: Output the inversion results.

[0090] During iterative training, the total loss function is continuously monitored. The trend of change. When the total loss function converges (i.e., the total loss value tends to stabilize and no longer decreases significantly), training stops. The criteria for convergence can be one or more of the following: reaching the preset maximum number of iterations; the decrease in the total loss value being less than a preset threshold for multiple consecutive steps (e.g., ...). The total loss value is less than the preset target value.

[0091] After training stops, the finally converged GU-CBAM-Unet inversion network model is forward-propagated once with the same input data to obtain the final output feature data, which is an Nz×Nx two-dimensional matrix. Each element in the matrix represents the logarithmic resistivity of the corresponding model grid cell. Taking the power of 10 of this matrix converts it into actual resistivity values ​​(in Ω·m). This two-dimensional resistivity matrix is ​​the high-resolution two-dimensional resistivity model of the target geological exploration area.

[0092] This two-dimensional resistivity model can be visualized as a color profile, with the horizontal axis representing the distance along the survey line and the vertical axis representing the depth. The color represents the magnitude of the resistivity value. Geophysicists can use this model to identify the location, shape, size, and resistivity characteristics of subsurface anomalies, providing quantitative information on the deep electrical structure for mineral resource exploration, geological structure research, and engineering geological surveys.

[0093] The implementation principle of this application's embodiments is as follows: The core principle of the proposed magnetotelluric two-dimensional inversion method lies in constructing a specifically improved GU-CBAM-Unet deep neural network and establishing an end-to-end differentiable mapping link between the network's output resistivity model and measured impedance data using a differentiable magnetotelluric physical forward modeling operator as a bridge. Based on this, a carefully designed total loss function, including data fitting terms and model smoothness constraints, drives the network to spontaneously search for a two-dimensional resistivity model that both fits the measured data well and possesses geological rationality, even with only single-sample, unlabeled measured data.

[0094] Replacing the batch normalization layer with a group normalization layer is a key technical guarantee for achieving unsupervised training with single samples, ensuring the stability of the feature normalization statistics under extremely small batch conditions. The introduction of the CBAM attention mechanism enhances the network's ability to extract and focus on key features of subsurface electrical anomalies. The differentiable forward modeling operator embeds the physical laws of geophysical forward problems into the network's training loop, providing a physical driving force for label-free training. The model roughness constraint term serves as a prior regularization method, effectively mitigating the ill-posedness of the inversion problem. The synergistic effect of these four technologies enables this method to achieve stable and high-resolution subsurface two-dimensional resistivity model inversion using only measured data obtained from a single exploration, without relying on labeled datasets.

[0095] The embodiments described in this specific implementation are preferred embodiments of this application and are not intended to limit the scope of protection of this application. Identical components are represented by the same reference numerals. Therefore, all equivalent changes made to the structure, shape, and principle of this application should be covered within the scope of protection of this application.

Claims

1. A two-dimensional magnetotelluric inversion method, characterized in that, Includes the following steps: Step S1: Obtain the measured magnetotelluric impedance data of the exploration area to be inverted, and normalize the measured impedance data to construct the network input data; at the same time, obtain the coordinates and frequency points of the resistivity model subdivision nodes of the exploration area. Step S2: Construct a GU-CBAM-Unet inversion network model; the GU-CBAM-Unet inversion network model adopts the Unet basic network structure, replaces the batch normalization layer in the basic network structure with the group normalization GN layer, and introduces the convolutional attention CBAM module at the downsampling path and skip connection; Step S3: Based on the differentiable magnetotelluric two-dimensional forward modeling operator, construct a total loss function that includes fitting error terms and model roughness constraint terms; Step S4: Input the network input data, resistivity model partitioning node coordinates and frequency point data as single samples into the GU-CBAM-Unet inversion network model, and drive the network to perform iterative training with the total loss function to optimize the weight parameters of the network model; Step S5: Stop training when the total loss function converges, and output the feature data of the finally converged GU-CBAM-Unet inversion network model as a high-resolution two-dimensional resistivity model of the target geological exploration area.

2. The magnetotelluric two-dimensional inversion method according to claim 1, characterized in that, In step S1, the magnetotelluric measured impedance data includes apparent resistivity data and impedance phase data in TM polarization mode and / or TE polarization mode.

3. The magnetotelluric two-dimensional inversion method according to claim 2, characterized in that, In step S1, the measured impedance data is normalized, specifically by using the following formula for the apparent resistivity data: Normalization is performed; for the impedance phase data, the formula is: Normalization is performed; among which, The normalized apparent resistivity To measure the apparent resistivity, and These are the maximum and minimum values ​​of the measured logarithmic resistivity, respectively. The normalized impedance phase, The measured impedance phase, and These represent the maximum and minimum values ​​of the measured impedance phase, respectively.

4. The magnetotelluric two-dimensional inversion method according to claim 3, characterized in that, When using TM polarization mode combined with TE polarization mode inversion, the steps to construct the network input data include: sequentially merging and splicing the normalized TM apparent resistivity, TM impedance phase, TE apparent resistivity and TE impedance phase to form a four-channel multidimensional tensor with a shape of 1×4×number of frequency points×number of observation points, which serves as the single-sample training input data for the network.

5. The magnetotelluric two-dimensional inversion method according to claim 1, characterized in that, In step S2, the group normalization GN layer does not rely on the batch dimension of the input data for statistics. Instead, it groups the channel features and calculates the mean and variance within each group for normalization to ensure the stability of the feature distribution of single-sample measured data during training iterations.

6. The magnetotelluric two-dimensional inversion method according to claim 1, characterized in that, In step S2, the Convolutional Attention (CBAM) module is a dual-path structure, consisting of a cascaded channel attention submodule and a spatial attention submodule. This module is used to suppress the geophysical background interference features extracted by the network and focus on the spatial geometric location features of deep electrical anomalies.

7. The magnetotelluric two-dimensional inversion method according to claim 6, characterized in that, The internal processing logic of the Convolutional Attention (CBAM) module is executed according to the following formula: F′=M c (F)⨂F(1) F′′=M s (F′)⨂F′(2) M c (F)= σ ( MLP ( AvgPool (F))+ MLP ( MaxPool (F)))(3) M s (F)= σ ( f 7×7 ([ AvgPool (F); MaxPool (F)]))(4) Wherein, the input feature map F∈R C×H×W F′ is the feature map after channel attention weighting, and M c ∈R C×1×1 For channel attention maps, M s ∈R 1×H×W F′′ is the spatial attention map, and F′′ is the final output feature map. σ It is the Sigmoid activation function. f 7×7 For 7×7 size convolution operations, MLP This represents a multilayer perceptron. AvgPool and MaxPool ∠ and ∠ represent average pooling and max pooling operations, respectively, and ⨂ represents element-wise multiplication.

8. The magnetotelluric two-dimensional inversion method according to claim 1, characterized in that, In step S3, the total loss function is calculated using the following formula: ;in, Total loss; This represents the fitting error term between the apparent resistivity forward response of the resistivity model output by the network and the measured normalized apparent resistivity. The term represents the fitting error between the impedance phase forward response of the resistivity model output by the network and the measured normalized impedance phase; m∈{TM,TE} represents the polarization mode. Θ represents the roughness constraint term of the model; Θ represents the two-dimensional resistivity model output by the network; ∇x and ∇z represent the gradients in the transverse and depth directions of the survey line, respectively; This indicates the calculation of the L2 norm; λ is the model smoothness weighting coefficient determined based on the environmental noise level.

9. A two-dimensional magnetotelluric inversion method according to claim 8, characterized in that, The fitting error term is calculated using the mean absolute error (MAE), and the specific formula is as follows: ; ;in, The differentiable magnetotelluric two-dimensional forward modeling operator is used to convert the two-dimensional resistivity model Θ output by the current iterative network into the corresponding theoretical impedance response in real time. and Indicates the measured magnetotelluric apparent resistivity and impedance phase; This indicates a normalization operation. This indicates the calculation of the mean absolute error.

Images