A method and system for magnetotelluric forward and inverse modeling based on neural operators

Abstract

This invention discloses a magnetotelluric (MT) forward and inverse modeling method and system based on neural operators, which falls within the field of geophysical exploration technology. It integrates the global feature extraction capability of Fourier neural operators with the continuous observation coordinate mapping mechanism of depth operator networks to construct a branch-backbone structured neural operator forward model. The branch network is used to extract the global medium feature vector of the subsurface two-dimensional resistivity model, while the backbone network uses observation location and frequency as continuous inputs to generate continuous coordinate basis functions. This enables the model to continuously predict the MT electromagnetic response at any observation location and frequency within the training range, reducing dependence on a fixed discrete output grid. Furthermore, this invention introduces a hybrid periodic implicit representation module composed of smooth activation functions and sinusoidal periodic activation functions into the backbone network to enhance the model's ability to characterize electromagnetic field attenuation and oscillation propagation features.

Description

Technical Field

[0001] This invention belongs to the field of magnetotelluric (MT) technology, specifically relating to a magnetotelluric (MT) forward and inverse method and system based on neural operators. Background Technology

[0002] Magnetotelluric (MT) is an important geophysical exploration method that uses natural alternating electromagnetic fields as a source to probe the electrical structure of the Earth's deep interior. It plays an irreplaceable role in mineral resource exploration and deep structural research. The interpretation of MT data is highly dependent on the reliability of the inversion results, while the accuracy and efficiency of the inversion are directly limited by the accuracy and computational speed of the forward modeling.

[0003] The core of MT forward modeling lies in solving the Helmholtz equations derived from Maxwell's equations. Traditional numerical simulation methods, such as finite difference, finite element, and finite volume methods, transform partial differential equations into large-scale linear equation systems by spatially discretizing the computational domain. These methods are based on physical laws and have high computational accuracy. However, their computational efficiency is highly dependent on the fineness of the mesh; finer meshes can improve forward modeling accuracy but significantly increase computational costs. This computational bottleneck is particularly prominent in nonlinear inversion problems that require hundreds or thousands of forward modeling calls.

[0004] To overcome the efficiency bottleneck of traditional numerical solutions, deep learning (DL) technology has been widely used to construct fast forward modeling surrogate models for medium-scale (MT) modeling. In recent years, advanced network architectures based on Transformer and attention mechanisms have further improved the feature extraction capabilities of these models; combined with mechanisms such as attention gating, DL has also shown great potential in 3D forward modeling. However, whether in 2D or 3D, the aforementioned conventional DL models essentially learn nonlinear mappings between finite-dimensional vector spaces, exhibiting severe "resolution dependence." These networks can only make predictions based on fixed observation parameters (such as measurement point locations and frequencies) used during training. Once the observation parameters change, training data must be regenerated and the network retrained, which severely limits their flexibility and applicability in practical exploration.

[0005] To overcome the "resolution dependence" of conventional deep learning models on position and frequency, neural operator techniques based on infinite-dimensional function space mappings have been introduced into the field of MT forward modeling. Neural operators aim to learn the continuous mapping patterns of underlying partial differential equations and have become a new generation of computational paradigms to replace traditional numerical solvers.

[0006] However, existing neural operator methods still face core challenges in MT forward modeling. There is an urgent need in this field for an efficient forward modeling surrogate model capable of continuous prediction at arbitrary observation locations and frequencies, in order to improve the accuracy, efficiency, and generalization ability of MT forward simulation and inversion calculations. Summary of the Invention

[0007] This invention aims to overcome the technical shortcomings of existing neural operator methods in 2D magnetotelluric (MT) forward modeling, significantly improve the model's accuracy in characterizing electromagnetic responses, break through the limitations of fixed grids, and enhance generalization performance under different measurement points and frequency sampling conditions. To this end, this invention provides a magnetotelluric forward and inverse modeling method and system based on neural operators. This method integrates the global feature extraction capability of the Fourier Neural Operator (FNO) with the continuous observation coordinate mapping mechanism of the DeepONet, enabling the neural operator forward modeling model to predict high-precision electromagnetic responses at any observation location and frequency within the training range. This breaks through the dependence of traditional deep learning models on fixed discrete grids, significantly improving the model's flexibility and applicability in practical exploration.

[0008] Specifically, this invention employs a branch-backbone neural operator forward model, where the branch network extracts global medium feature vectors from the subsurface two-dimensional resistivity model, and the backbone network takes observation location and frequency as continuous inputs, introducing a hybrid periodic neural representation dual-branch fusion mechanism. Thus, the neural operator forward model can continuously predict the MT electromagnetic response at any observation location and frequency within the training range, and can be embedded as a differentiable forward modeling engine into the inversion framework to achieve rapid forward and inversion of magnetotelluric data.

[0009] Based on the above technical approach, this invention further optimizes the technology by introducing a hybrid smoothing-periodic activation mechanism into the MT operator learning process. This involves combining the SiLU smoothing activation function with a sinusoidal periodic activation function to simultaneously characterize the smooth decay trend and periodic oscillation characteristics in the MT response. Specifically, the SiLU component characterizes the smooth decay trend of the MT response with variations in spatial location and frequency, while the sinusoidal periodic component characterizes the periodic oscillation characteristics of the electromagnetic field phase with variations in frequency and propagation path. Ultimately, this leads to the preferred model of this invention: the Extended Fourier Hybrid Periodic Implicit Neural Operator (EFHPN).

[0010] Therefore, the present invention provides the following technical solution:

[0011] A magnetotelluric forward modeling method based on neural operators includes the following steps:

[0012] Step 1: Obtain the underground two-dimensional resistivity model ρ(y,z) and the horizontal coordinate y and frequency f of the observation point;

[0013] Step 2: Input the data obtained in Step 1 into the neural operator forward model. The neural operator forward model adopts the branch-backbone dual network architecture of the DeepONet deep operator network, which includes a backbone network, branch networks, and an operator fusion module. The Fourier neural operator FNO is introduced into the branch network. The branch network is used to extract the global medium feature vector b from the two-dimensional resistivity model ρ(y,z) of the subsurface. k k is the feature channel index; the backbone network introduces a hybrid periodic neural representation dual-branch fusion mechanism, and the backbone network is used to generate continuous coordinate basis functions t based on the continuous observation coordinates (y,f). k ;

[0014] The operator fusion module will integrate the continuous coordinate basis function t k With global medium eigenvector b k By fusing tensor inner products, the electromagnetic response at the corresponding observation location and frequency is output, thus obtaining the forward modeling result;

[0015] The electromagnetic response includes TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase.

[0016] It should be understood that the Fourier Neural Operator (FNO), with its efficient frequency domain feature extraction capabilities, has become the most representative basic architecture. By extending the FNO architecture, prediction of arbitrary measurement point locations is achieved; further integration with the DeepONet mechanism enhances the generalization ability to untrained frequencies. Thanks to the continuous mathematical framework of operator learning, these models possess strong generalization capabilities, enabling them to quickly predict electromagnetic responses at any given spatial location and frequency.

[0017] Optionally, the branch network has a sequentially connected up-dimensional layer P, a number of cascaded Fourier neural operator layers, and a down-dimensional layer Q;

[0018] Each Fourier neural operator layer employs the Fourier neural operator FNO, comprising two parallel branches. The first branch follows the sequence: "Fourier transform F → frequency domain spectral weighting modulation R → inverse Fourier transform F". -1 The frequency domain convolution paradigm is used; the second branch uses a local linear transformation W as a bypass connection; finally, the outputs of the two branches are summed and fused to obtain the latent feature tensor.

[0019] The latent feature tensor obtained from the current layer is used as the input latent feature tensor of the next Fourier neural operator layer, and the high-dimensional latent feature tensor output from the last Fourier neural operator layer is then input to the dimensionality reduction layer Q.

[0020] Optionally, the backbone network includes a normalization layer, multiple hidden fully connected layers, and an output layer connected in sequence.

[0021] The normalization layer is used to normalize the horizontal coordinate y and frequency f of the observation point;

[0022] The multiple hidden fully connected layers are used to perform high-dimensional nonlinear feature mapping on the normalized continuous observation coordinates (y,f), and the hidden fully connected layers adopt a hybrid periodic neural representation dual-branch fusion mechanism, that is, one branch adopts a periodic activation function and the other branch adopts a nonlinear activation function.

[0023] The output layer is used to generate continuous coordinate basis functions t. k .

[0024] Optionally, the hybrid periodic neural network represents a bibranch fusion mechanism in which the two branches are a smooth branch and a sinusoidal periodic branch, respectively;

[0025] The smooth branch uses the SiLU activation function, and the sinusoidal periodic branch uses the sinusoidal activation function. The outputs of the two branches are summed and fused element-wise. The corresponding hybrid periodic neural representation is as follows:

[0026] ;

[0027] in, This is a branch of the SiLU activation function, used to characterize the smooth decay trend and background changes in the electromagnetic response; W1 and W2 are the sinusoidal activation function branches used to characterize the periodic oscillations and phase changes in the electromagnetic response; W1 and W2 are learnable weight matrices, and ω0 is a hyperparameter controlling the frequency scale of the sinusoidal function. denoted as the hybrid activation function in the hidden fully connected layer of the backbone network; x represents the input of the hidden fully connected layer.

[0028] It should be understood that when an electromagnetic field propagates in an underground conductive medium, it exhibits both periodic phase changes and energy attenuation with increasing depth. Traditional neural operators using conventional activation functions such as ReLU, Tanh, and GELU are easily affected by spectral bias, making it difficult to accurately characterize both the rapidly changing phase response and the slowly changing amplitude attenuation simultaneously. The technical solution of this invention, through the aforementioned design, can simultaneously characterize the periodic oscillation characteristics and energy attenuation law in the MT electromagnetic response.

[0029] This invention also provides a magnetotelluric inversion method based on neural operators, comprising the following steps:

[0030] Construct an inversion framework based on a neural operator forward model to accelerate the inversion process: use the pre-trained neural operator forward model as the forward modeling engine;

[0031] Define the objective function and the initial model parameters m0;

[0032] Acquisition and preprocessing of observation data: Observation data d of the inversion area were acquired using MT detection equipment. obs ;

[0033] Iterative inversion: based on initial model parameters m0, a neural operator forward model as the forward modeling engine, and observation data d. obs Perform iterative inversion and update the model parameters m;

[0034] Here, the model parameter m is the target parameter model of the inversion area, which represents the complex underground electrical structure; the model parameter m obtained at the end of the iteration is used as the inversion result of the underground electrical structure of the inversion area.

[0035] Optionally, the inversion process employs L-BFGS regularized inversion, and the objective function is constructed using the Tikhonov regularization framework. The L-BFGS regularized inversion process is as follows:

[0036] Based on the initial model parameters m0, the neural operator forward model is invoked to calculate the model parameters m corresponding to the k-th iteration. k The corresponding predicted response;

[0037] Based on predicted response and observation data d obs Construct the objective function And use automatic differentiation to calculate the objective function. Relative to m k gradient k is the iteration round marker;

[0038] Gradient based on the objective function The L-BFGS algorithm is used to update the search direction p. k The step size α is determined by line search. k This updates the model parameters corresponding to the (k+1)th iteration. Enter the iterative operation;

[0039] The iteration terminates when the gradient norm is less than a preset threshold or the change in the objective function is less than a preset threshold.

[0040] The model parameters m at the end of the iteration inv Output as the target parameters of the inversion region.

[0041] The present invention also provides a forward modeling system based on the above-described forward modeling method, comprising:

[0042] The model building module is used to build a neural operator forward model that includes a branch network, a backbone network, and an operator fusion module.

[0043] The data acquisition module is used to acquire the underground two-dimensional resistivity model ρ(y,z) and the horizontal coordinate y and frequency f of the observation point;

[0044] The forward modeling module is used to input the subsurface two-dimensional resistivity model ρ(y,z) into the branch network of the neural operator forward model, and input the horizontal coordinate y and frequency f of the observation point into the backbone network of the neural operator forward model. The operator fusion module realizes the continuous mapping from the subsurface resistivity distribution to the surface electromagnetic response.

[0045] The electromagnetic response includes TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase.

[0046] This invention also provides an inversion system, comprising:

[0047] The inversion model building module is used to use the pre-trained neural operator forward model as the forward modeling engine;

[0048] The observation module is used to collect observation data of the inversion area through the MT detection equipment. obs ;

[0049] The inversion module is used to perform forward modeling based on initial model parameters m0, a neural operator forward model serving as the forward modeling engine, and observation data d. obs The model parameters m are updated through iterative inversion to finally obtain the complex underground electrical structure of the inverted area.

[0050] The present invention also provides a computer terminal, comprising: one or more processors and a memory storing a computer program;

[0051] The processor invokes a computer program to implement the steps of the above-mentioned neural operator-based magnetotelluric forward modeling method or the steps of the above-mentioned neural operator-based magnetotelluric inversion method.

[0052] The present invention also provides a computer-readable storage medium storing a computer program, which is invoked by a processor to implement: the steps of the above-described magnetotelluric forward modeling method based on neural operators, or the steps of the above-described magnetotelluric inversion method based on neural operators.

[0053] Compared with the prior art, the present invention achieves the following progress and effects:

[0054] (1) Achieving continuous prediction of location and frequency, breaking through grid dependence. This invention, through the continuous observation coordinate input mechanism of DeepONet, enables the forward model to predict high-precision electromagnetic responses at any observation location and frequency within the training range, breaking through the dependence of traditional deep learning models on fixed discrete grids, and significantly improving the flexibility and applicability of the model in actual exploration; in particular, the technical solution of this invention integrates the global feature extraction capability of Fourier Neural Operator (FNO) and the continuous observation coordinate mapping mechanism of DeepONet. The FNO branch network is responsible for extracting the global medium feature vector of the underground two-dimensional resistivity model, and the DeepONet backbone network is responsible for establishing the continuous observation coordinate mapping from observation location and frequency to electromagnetic response. After the two are fused by the operator, the model has both global structure perception capability and continuous prediction capability of arbitrary coordinates, thereby breaking through the fixed grid limitation and improving the generalization performance under different measurement point and frequency sampling conditions.

[0055] (2) A hybrid smooth-periodic activation mechanism is introduced to enhance the characterization capability of electromagnetic response. In particular, the improved technical solution of this invention introduces a dual-branch fusion mechanism of HPN composed of the SiLU activation function and the sinusoidal periodic activation function in the backbone network. The SiLU branch is used to characterize the smooth decay trend and background changes in the MT response, while the sinusoidal periodic branch is used to characterize the periodic oscillations and phase change characteristics in the electromagnetic response. After fusion, the two can simultaneously characterize the energy decay law and phase periodic changes in the electromagnetic field propagation process, thereby alleviating the spectral bias problem of conventional activation functions and improving the accuracy and stability of forward modeling prediction;

[0056] (3) Construct a differentiable forward surrogate model and integrate it into the L-BFGS inversion framework. This invention utilizes automatic differentiation technology to efficiently calculate the gradient of the objective function, enabling the forward model, especially EFHPN, to be seamlessly embedded as a differentiable module into gradient-based optimization algorithms. This achieves fast and high-precision inversion of complex two-dimensional resistivity structures, significantly reducing the computational cost of traditional numerical inversion methods. Attached Figure Description

[0057] Figure 1 This is a schematic diagram of electromagnetic signal propagation in underground media, where, Figure 1 (a) is a schematic diagram of the attenuation and oscillation of an electromagnetic field propagating in an underground medium, and Figure 1 (b) is a schematic diagram of the characterization of the smooth component, sinusoidal component and mixed activation component;

[0058] Figure 2 A schematic diagram of the MT fast forward modeling network architecture for the extended Fourier hybrid periodic implicit neural operator EFHPN.

[0059] Figure 3A schematic diagram of the magnetotelluric two-dimensional resistivity model generated by forward mesh partitioning and Gaussian random field generation, wherein (a) is a schematic diagram of the mesh partitioning of the MT computational domain, with black lines representing the mesh and black boxes representing the core region; (b) is the resistivity model of the core region.

[0060] Figure 4 The diagram shows the changes in training loss and test loss during the training process of four different forward modeling networks. (a) is a schematic diagram of test loss; (b) is a schematic diagram of training loss; and (c) is a schematic diagram of changes in training learning rate.

[0061] Figure 5 The figures show a comparison of the forward prediction results of four neural operator forward models on the test sample, where (a) is the resistivity model of the test sample; (b1), (b2), (b3), and (b4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the FDM benchmark solution, respectively; (c1), (c2), (c3), and (c4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the EFNO prediction results, respectively; and (d1), (d2), and (d3) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the EFNO prediction results, respectively. (d1), (e2), (e3), (e4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase corresponding to the EFPN prediction results, respectively; (e1), (e2), (e3), (e4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase corresponding to the EFKAN prediction results, respectively; (f1), (f2), (f3), (f4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase corresponding to the EFHPN prediction results, respectively.

[0062] Figure 6 The figures show a comparison of the apparent resistivity and phase prediction errors of four neural operator forward models on the test samples. (a1), (a2), (a3), and (a4) represent the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the EFNO prediction error, respectively; (b1), (b2), (b3), and (b4) represent the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the EFPN prediction error, respectively; (c1), (c2), (c3), and (c4) represent the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the EFKAN prediction error, respectively; and (d1), (d2), (d3), and (d4) represent the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase in the EFHPN prediction error, respectively.

[0063] Figure 7The figure shows a comparison of the forward response prediction results of four neural operator forward models at different observation locations. In the figure, (a1), (a2), (a3), and (a4) are the apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase of the forward response at a location of -15km, respectively; and (b1), (b2), (b3), and (b4) are the apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase of the forward response at a location of 15km, respectively.

[0064] Figure 8 The figure shows a comparison of the forward response prediction results of four neural operator forward models at different frequencies. Among them, (a1), (a2), (a3), and (a4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase corresponding to the forward response at 0.02 Hz, respectively; (b1), (b2), (b3), and (b4) are the TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase corresponding to the forward response at 500 Hz, respectively.

[0065] Figure 9 Comparison of forward modeling acceleration L-BFGS inversion results for different neural operators: (a) results of theoretical model, (b) inversion results of EFNO+L-BFGS, (c) inversion results of EFPN+L-BFGS, (d) inversion results of EFKAN+L-BFGS, and (e) inversion results of EFHPN+L-BFGS. Detailed Implementation

[0066] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of the invention. The technical features involved in the various embodiments of the invention described below can be combined with each other as long as they do not conflict with each other.

[0067] It should be noted that although functional modules are divided in the device schematic diagram and a logical order is shown in the flowchart, in some cases, the steps shown or described may be performed in a different order than those in the device or the flowchart. The terms "first," "second," etc., used in the specification, claims, and drawings are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. Unless otherwise defined, all technical and scientific terms used herein have the meanings commonly understood by one of ordinary skill in the art to which this application pertains. The terminology used herein is only for describing embodiments of this application and is not intended to limit this application.

[0068] This invention proposes a magnetotelluric (MT) forward and inverse modeling method and system based on neural operators, constructing a fast forward modeling model that combines prediction accuracy, continuous generalization ability, and differentiability. This method integrates the global feature extraction capability of the Fourier Neural Operator (FNO) with the continuous observation coordinate mapping mechanism of the DeepONet operator network, forming a branch-backbone structured neural operator forward model. Specifically, the branch network extracts the global medium feature vector of the subsurface two-dimensional resistivity model, while the backbone network generates continuous coordinate basis functions based on the observation location and frequency, thereby enabling prediction of the MT electromagnetic response at any observation location and frequency within the training range.

[0069] This invention also proposes introducing a hybrid periodic neural representation (HPN) composed of a SiLU smooth activation function and a sinusoidal periodic activation function into the backbone network. The SiLU branch characterizes the smooth decay trend and background changes in the electromagnetic response, while the sinusoidal periodic branch characterizes the periodic oscillations and phase changes in the electromagnetic response. Through this structure, this invention enhances the expressive power of the forward model for complex MT responses and allows the pre-trained forward model to be embedded as a differentiable forward modeling engine within the L-BFGS regularized inversion framework, enabling rapid calculation of the inversion gradient using automatic differentiation.

[0070] When this invention is applied to two-dimensional magnetotelluric inversion, the inversion efficiency is better than that of traditional nonlinear conjugate gradient inversion, and the inversion results are highly consistent with the shallow seismic interpretation profile, thus verifying the effectiveness and feasibility of this invention.

[0071] Magnetotelluric two-dimensional forward modeling problem

[0072] Regarding the two-dimensional forward modeling of magnetotellurics (MT), to facilitate understanding of this invention, a brief explanation of the two-dimensional MT forward modeling problem is provided below. The theoretical basis of MT is the propagation and attenuation law of electromagnetic fields in subsurface media. Assume the time harmonic factor is e. -iωt Under the conditions of quasi-static approximation and passive region, the frequency domain Maxwell's equations can be expressed as:

[0073] ;

[0074] In the formula, E is the electric field strength (V / m), H is the magnetic field strength (A / m), and μ is the permeability of the medium (usually taken as the permeability of free space). H / m), σ is the electrical conductivity of the underground medium (S / m), the resistivity of the underground medium ρ=1 / σ(Ω.m), ω=2πf is the angular frequency, and f is the exploration frequency.

[0075] For a typical two-dimensional geoelectric structure, assuming the geological structure extends infinitely along the strike (x-axis), the partial derivatives of all physical quantities in the strike direction are zero. Furthermore, the dielectric conductivity is solely a function of the cross-sectional coordinates (y,z), σ(y,z). Under these conditions, the electromagnetic field can be decomposed into two independent polarization modes: the transverse electric mode TE and the transverse magnetic mode TM.

[0076] In TE mode, the electric field retains only the directional component E. x The magnetic field includes a cross-sectional component H. y and H z The partial differential control equations for the TE mode can be derived from Maxwell's equations:

[0077] ;

[0078] In TM mode, the magnetic field retains only the directional component H. x The electric field includes a cross-sectional component E. y and E z Similarly, the partial differential control equations for the TM mode can be derived:

[0079] ;

[0080] The geological exploration target area is divided into a core target area and a peripheral extension area. For the core target area containing geological anomalies, in some implementations, a uniform grid is used horizontally, while a non-uniform grid with gradually increasing grid size is used vertically to accommodate the skin effect of electromagnetic field energy attenuation with depth. In some implementations, a non-uniform grid with a stretching factor greater than 1 can be used for the peripheral extension area. To balance computational accuracy and efficiency, the finite difference method can be used to solve the aforementioned governing equations, thereby extracting the surface electromagnetic response and calculating the apparent resistivity and phase.

[0081] After acquiring the electromagnetic field components at the Earth's surface, the surface impedance component Z in the TE and TM modes... TE and Z TM They are defined as follows:

[0082] ;

[0083] Based on the above impedance, the apparent resistivity (ρ) can be further calculated. a ) and phase (φ):

[0084] ;

[0085] In the formula, , Apparent resistivity in TE mode and TM mode, respectively; , , respectively, represent the impedance phase in TE mode and TM mode; Im and Re represent the imaginary and real parts of the complex number, respectively.

[0086] Therefore, the two-dimensional MT forward modeling problem can be formalized as an operator mapping from the subsurface two-dimensional resistivity model ρ(y,z) to the surface electromagnetic response. :

[0087] ;

[0088] In the formula, ρ(y,z) represents the subsurface two-dimensional resistivity model, (y,f) represents the horizontal position and frequency of the observation point, and the output of the two-dimensional MT forward modeling problem is the TE mode apparent resistivity at the corresponding observation position and frequency. TE mode phase TM mode apparent resistivity Phase with TM mode .

[0089] This invention does not directly solve the above partial differential equations, but instead trains a neural operator to approximate the operator mapping. In the following embodiments, the Extended Fourier Mixed Periodic Implicit Neural Operator (EFHPN) is used as the best example for illustration. During the training phase, a large number of "two-dimensional resistivity model-electromagnetic response" data pairs are generated using the finite difference method; during the prediction phase, the neural operator forward model directly outputs the electromagnetic response through forward propagation, thereby significantly reducing the computation time of forward modeling.

[0090] Neural operator forward model:

[0091] Conventional deep learning networks can typically only learn mappings between finite-dimensional discrete vectors, and their generalization ability is limited by the grid resolution during training. The core advantage of neural operators lies in their ability to learn nonlinear mappings between infinite-dimensional function spaces, thus breaking through the limitations of fixed discrete grids and becoming a type of computational model to replace traditional numerical solvers for partial differential equations. This invention introduces two representative neural operator architectures: DeepONet (deep operator network) and the Fourier neural operator FNO (Fourier No. 1).

[0092] DeepONet, a deep operator network, employs a branch-backbone dual network architecture to decouple the input function from the output coordinates. Its operator expression can be written as:

[0093] ;

[0094] Where ρ is the input two-dimensional resistivity model, i.e., ρ(y,z) as described above. To represent a branch network, b k Let r1, r2, r be the k-th global dielectric feature vector extracted by the branch network from the two-dimensional resistivity model ρ. m Let r be the 1st, 2nd, and mth discrete spatial sampling points or grid nodes in the two-dimensional resistivity model ρ, and r i =(y i ,zi ), y i For sampling points or grid nodes r i horizontal coordinate, z i For sampling points or grid nodes r i The depth coordinates; ρ(r1), ρ(r2), ..., ρ(r m ) represents the two-dimensional resistivity model ρ at sampling points or grid nodes r1, r2, ..., r m The resistivity value at that location; is the k-th continuous coordinate basis function generated by the backbone network based on the horizontal position y and frequency f of the observation point; p is the latent feature dimension, which is determined by the dimension of the final output layer of the branch network and the backbone network.

[0095] Because this invention employs a fusion method that multiplies and sums the global medium feature vector of the branch network with the continuous coordinate basis functions of the backbone network term by term, the number of features output by the branch network remains consistent with the number of continuous coordinate basis functions output by the backbone network. This architecture enables the network to make predictions at arbitrary horizontal measurement point locations and frequencies within a preset domain, thereby reducing dependence on a fixed discrete output grid.

[0096] Based on the mathematical model described above, the neural operator forward model comprises a backbone network, branch networks, and an operator fusion module. The branch networks employ a Fourier neural operator (FNO) architecture to extract global medium feature vectors from the two-dimensional discrete resistivity model ρ(y,z). The backbone network uses a hybrid periodic neural representation architecture to generate continuous coordinate basis functions from the continuous observation coordinates (y,f). The operator fusion module fuses the two in the latent space to output the corresponding MT electromagnetic response.

[0097] Branch network structure:

[0098] The branch network consists of a sequentially connected up-dimensional layer P, several cascaded Fourier neural operator layers, and a down-dimensional layer Q.

[0099] The dimension-upgrading layer P maps the input single-channel two-dimensional resistivity model ρ(y,z) to a high-dimensional latent feature space, obtaining the initial latent feature tensor. Each Fourier neural operator layer includes two parallel branches: the first branch follows the pattern "Fourier transform F → frequency domain spectral weight modulation R → inverse Fourier transform F". -1 The first branch uses a frequency-domain convolutional paradigm to extract global features from the latent feature tensor in the frequency domain. The second branch uses a local linear transformation W as a bypass connection to preserve local spatial information and anomalous body boundary details in the input latent feature tensor. The outputs of the two branches are summed and fused, and then passed through a nonlinear activation function to obtain the next latent feature tensor, which is then input into the next Fourier neural operator layer. The output of the last Fourier neural operator layer is a high-dimensional latent feature tensor, which is then input into the dimensionality reduction layer Q.

[0100] In this embodiment, the initial latent feature tensor after the dimensionality-upgrading layer P is iteratively updated and its features extracted through six layers of Fourier neural operators. During the frequency domain convolution, the branch network truncates the Fourier modes and retains the first 18 low-frequency modes to extract the main global structural information in the subsurface two-dimensional resistivity model and reduce computational complexity. Finally, the dimensionality-reduction layer Q compresses the high-dimensional latent features into a global medium feature vector b. k This is used for operator fusion with the coordinate-related basis functions (continuous coordinate basis functions) output by the backbone network, thereby predicting the MT electromagnetic response (TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase). It should be understood that the aforementioned feature dimensions, number of Fourier neural operator layers, and number of preserved modes are exemplary values ​​for this embodiment and do not constitute a limitation on the scope of protection of this invention. The dimensionality-upgrading layer P and the dimensionality-reducing layer Q can be selected from conventional network layers used for dimensionality upgrading and dimensionality reduction.

[0101] Fourier neural operator layer:

[0102] The Fourier neural operator (FNO) utilizes the convolution theorem to transform the global integral operator into frequency-domain spectral weight modulation in the Fourier domain, thereby efficiently extracting global spatial features. The inter-layer update rules for each Fourier neural operator layer are as follows:

[0103] ;

[0104] In the formula, F and F -1 σ represents the Fast Fourier Transform and its inverse transform, respectively; σ is the nonlinear activation function; R is the learnable spectral weight transformation in the Fourier domain; W is the local linear transformation, i.e., the corresponding convolution operator, which serves as a bypass connection to preserve the detailed information of the anomalous body boundary. For the input subsurface two-dimensional resistivity model, v j Let be the input latent feature tensor of the j-th Fourier neural operator layer, representing the resistivity model. The high-dimensional feature representation in the j-th Fourier neural operator layer, v j+1 The latent feature tensor output by the j-th Fourier neural operator layer is used as the input latent feature tensor of the next (j+1)-th Fourier neural operator layer. By stacking multiple layers of FNO, the branch network can capture long-range spatial dependencies with logarithmic parameter complexity.

[0105] Backbone network structure:

[0106] The backbone network takes as input the set observation data coordinates and frequencies (y, f), and outputs the continuous coordinate basis function t. k The backbone network consists of a normalization layer, multiple hidden fully connected layers, and an output layer connected in sequence.

[0107] The normalization layer is used to normalize the horizontal coordinate y and frequency f of the observation point. In this embodiment, the frequency f is preferably first logarithmically transformed and then normalized to the [-1,1] interval to reduce the impact of different physical quantities on the stability of network training.

[0108] Multiple hidden fully connected layers are used to perform high-dimensional nonlinear feature mapping on the normalized continuous observation coordinates (y, f). In this embodiment, the backbone network can be a 4-layer fully connected network, with the number of neurons in each layer being (2→128→256→512→1024), where the two dimensions of the input correspond to the horizontal position y and frequency f of the observation point. Each hidden fully connected layer adopts a two-branch fusion mechanism, that is, it has two branches, each using a different activation function. The outputs of the two branches are summed element-wise and fused, then passed through LayerNorm normalization and PReLU activation, and used as the input to the next hidden fully connected layer. The output of the last hidden fully connected layer is input to the output layer, and the output layer obtains the continuous coordinate basis function t. k .

[0109] The preferred embodiment of this invention employs a hybrid periodic neural representation (HPN) dual-branch fusion mechanism in the hidden fully connected layer, where the smooth branch uses the SiLU activation function and the sinusoidal periodic branch uses the sinusoidal activation function, as will be described in detail below. In other embodiments, the activation functions described above can be replaced with other nonlinear or periodic activation functions as needed, but the combination of the SiLU activation function and the sinusoidal activation function is the preferred option.

[0110] Operator fusion module: The core link realizing the mapping from the physical domain (resistivity model) to the response domain (electromagnetic response). The global medium feature vector output by the branch network and the continuous coordinate basis functions generated by the backbone network are dually fused in the latent space through tensor inner product. This design assigns independent feature channels to each electromagnetic response component to learn the polarization differences between TE and TM modes, while sharing the same backbone network to ensure the consistency of coordinate encoding. The entire forward modeling process is fully differentiable end-to-end, enabling it to be seamlessly embedded as an independent forward modeling surrogate model into the gradient-based L-BFGS inversion optimization framework.

[0111] Hybrid periodic neural networks (HPNs) are neural networks that combine a smooth activation function and a periodic activation function in a two-branch fusion mechanism.

[0112] In some embodiments, taking into full account that existing networks generally use conventional activation functions such as ReLU and Tanh, whose mathematical expressions are difficult to simultaneously fit the smooth decay trend and periodic oscillation characteristics in electromagnetic response, resulting in the model's limited ability to characterize complex MT responses, this invention further proposes a hybrid periodic neural representation dual-branch fusion mechanism, where one branch is a nonlinear activation function, preferably a smooth activation function, and the other branch is a periodic activation function.

[0113] Although DeepONet and FNO possess strong operator mapping capabilities, conventional activation functions often suffer from spectral bias, meaning the network tends to preferentially learn low-frequency smooth components while lacking the ability to represent rapidly changing or periodically oscillating components. When magnetotelluric responses propagate in subsurface media, they typically exhibit both periodic phase changes and amplitude attenuation characteristics, making it difficult for conventional activation functions to naturally fit this physical law (e.g., ...). Figure 1 (as shown in a).

[0114] To simultaneously capture the smooth background trend and periodic oscillation characteristics of the electromagnetic response, this invention also proposes a hybrid periodic neural representation (HPN) dual-branch fusion mechanism, which combines a smooth activation function with a sinusoidal periodic activation function. The feature mapping of the hybrid periodic neural representation (HPN) dual-branch fusion mechanism is defined as follows:

[0115] ;

[0116] Among them, the first item This is a branch of the SiLU activation function, used to capture the smooth decay trend and background changes in the electromagnetic response; the second term A sinusoidal activation function branch is used to characterize the periodic oscillations and phase changes in the electromagnetic response; W1 and W2 are learnable weight matrices; ω0 is a hyperparameter controlling the frequency scale of the sinusoidal function. In this embodiment, the input layer of the backbone network is set to ω0=45, and the subsequent hidden fully connected layers are set to ω0=1. denoted as the hybrid activation function in the hidden fully connected layer of the backbone network, used to enhance the network's ability to jointly express smooth decay features and periodic oscillation features, and x represents the input of the hidden fully connected layer.

[0117] In the preferred embodiment of this invention, the backbone network employs a Hybrid Periodic Neural Representation (HPN) architecture, specifically designed to handle continuous observation coordinate inputs (y, f). This design reduces the "resolution dependency" of conventional deep learning models, enabling the forward surrogate model to make continuous predictions at any location and frequency within the training range. This mechanism combines a smooth SiLU activation function with a sinusoidal activation function; the former is used to fit the decay trend and low-frequency background field in the electromagnetic response, while the latter is used to characterize the periodic oscillations and phase change features in the electromagnetic response (e.g., ...). Figure 1 (as shown in b).

[0118] EFHPN overall structure:

[0119] Based on the aforementioned network components, the preferred embodiment of this invention employs the Extended Fourier Hybrid Periodic Implicit Neural Operator (EFHPN). This EFHPN integrates the DeepONet deep operator network and the Fourier Neural Operator (FNO). Specifically, it uses the branch-backbone dual network architecture of DeepONet as the main model architecture, with the Fourier Neural Operator (FNO) incorporated into the branch networks. The hidden fully connected layers of the backbone network employ the Hybrid Periodic Neural Representation (HPN) architecture. This integrates the frequency domain global feature extraction capability of the Fourier Neural Operator (FNO), the continuous observation coordinate mapping mechanism of the DeepONet deep operator network, and the hybrid periodic representation capability of the Hybrid Periodic Neural Representation (HPN) (e.g.,...). Figure 2 (As shown). The overall mapping of EFHPN can be represented as:

[0120] ;

[0121] Where m and n represent the number of discrete grids or sampling points in the horizontal and depth directions of the input two-dimensional resistivity model, respectively. k The global dielectric feature vector extracted by the branch network FNO from the discrete resistivity model; t k The backbone network generates continuous coordinate basis functions based on the normalized continuous observation coordinates (y,f); the two are fused in the latent space through tensor inner product to obtain four electromagnetic response components at the corresponding observation positions and frequencies: TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase.

[0122] It should be understood that EFHPN is a preferred embodiment of the present invention. In other alternative embodiments, the sinusoidal periodic branches in the hidden fully connected layers of the backbone network can be replaced with cosine functions, sine-cosine combination functions, or other activation functions with periodic expressive capabilities; the SiLU activation function in the smoothing branch can also be replaced with nonlinear activation functions such as Tanh, GELU, Swish, ELU, or ReLU. The combination of SiLU and sinusoidal activation functions is a preferred scheme, as it can better characterize both the smooth decay characteristics and periodic oscillation characteristics in the MT response simultaneously. It also maintains the differentiability of the forward modeling network, facilitating seamless integration into the L-BFGS inversion framework. It should also be understood that forward models without incorporating periodic sinusoidal activation functions are still feasible techniques of the present invention and fall within the protection scope of the present invention.

[0123] Based on the above theoretical explanation, this invention provides a magnetotelluric forward modeling method based on neural operators, comprising the following steps:

[0124] Step 1: Obtain the subsurface two-dimensional resistivity model ρ(y,z) and the horizontal coordinates y and frequency f of the observation points. It should be understood that in practical applications, the forward modeling prediction range can be set according to the actual geological exploration needs.

[0125] Step 2: Input the data obtained in Step 1 into the neural operator forward model. The neural operator forward model adopts the branch-backbone dual network architecture of DeepONet, which includes a backbone network, branch networks and operator fusion module.

[0126] In this process, the Fourier neural operator FNO is introduced into the branch network. The two-dimensional resistivity model ρ(y,z) of the subsurface is input into the branch network to output the global medium feature vector b. k k is the feature channel index;

[0127] The backbone network incorporates a hybrid periodic neural representation (HPN) architecture, where one branch uses a periodic activation function and the other uses a nonlinear activation function. The horizontal coordinate y and frequency f of the observation point are input into the backbone network to generate a continuous coordinate basis function t. k ;

[0128] The operator fusion module will integrate the global medium feature vector b k With continuous coordinate basis function t k The electromagnetic response is output by fusing through tensor inner product within the hidden space. This electromagnetic response includes the apparent resistivity of TE mode, the phase of TE mode, the apparent resistivity of TM mode, and the phase of TM mode.

[0129] EFHPN accelerates L-BFGS regularized inversion:

[0130] This invention can also apply the above-mentioned neural operator forward model to MT inversion. The following will use the extended Fourier hybrid periodic implicit neural operator EFHPN applied to L-BFGS regularized inversion as an example. It should be understood that any neural operator forward model that does not deviate from the concept of this invention can be applied to MT inversion with reference to the following ideas.

[0131] MT inversion is essentially a nonlinear, ill-conditioned optimization problem. To obtain a stable subsurface resistivity distribution that conforms to geological laws, this embodiment employs a Tikhonov regularization framework to construct the objective function. :

[0132] ;

[0133] in, For data fitting terms, Here, d represents the model regularization term, m is the model parameter for the region to be inverted, such as resistivity model parameters; obs For MT observation data; F EFHPN(m) represents the response data predicted by the EFHPN forward model; W d A weighted matrix for the data; W m m is the model regularization matrix; ref λ is the reference model; λ is the regularization parameter used to balance the data fitting term and the model smoothing term.

[0134] To minimize the objective function, the finite-memory quasi-Newton method L-BFGS iterative update process can be expressed as:

[0135] ;

[0136] Where k is the iteration number, α k To satisfy the search step size of the Wolfe criterion, p k This determines the search direction for model updates.

[0137] The L-BFGS algorithm does not directly compute and store the complete Hessian matrix, but instead utilizes the parameter difference s from the most recent finite number of iterations. k =m k+1 -m k and gradient difference y k =g k+1 -g k (in The search direction is obtained by implicitly approximating the product of the Hessian matrix and the gradient using a double-loop recursive algorithm. This mechanism enables L-BFGS to achieve faster convergence speed with lower memory usage. Wherein, m k+1 ,m k The resistivity model parameters m and g are updated for the (k+1)th and kth iterations. k+1 ,g k Let be the gradient of the objective function in the (k+1)th and kth iterations.

[0138] However, regardless of the gradient-based inversion algorithm used, the core computational bottleneck lies in the gradient g of the objective function in each iteration. k To find g k The expression is:

[0139] ;

[0140] in, Let g be the sensitivity matrix (Jacobi matrix). In the traditional numerical solution framework, the gradient g is calculated. k Typically, it is necessary to solve the adjoint forward modeling problem to obtain the results. The product of the residual vector and the matrix is ​​computationally very expensive, and T is the matrix transpose symbol.

[0141] This invention significantly improves the efficiency of forward modeling by introducing a neural operator forward model, while the Jacobian matrix can be quickly obtained through automatic differentiation techniques within a deep learning framework. In particular, this embodiment introduces a fully differentiable EFHPN fast forward model, which can complete the forward modeling solution in milliseconds, and the Jacobian matrix can be quickly obtained through automatic differentiation techniques within a deep learning framework. Specifically, given any vector v, The calculation can be completed with only one backpropagation, without the need for explicit storage or solving of the sensitivity matrix. This strategy can effectively avoid the high computational cost of forward modeling and gradient solving in traditional L-BFGS inversion, achieving truly fast and high-precision inversion.

[0142] In summary, the L-BFGS regularized inversion method replaces the traditional finite difference forward solver with a pre-trained neural operator forward model. For example, the EFHPN-accelerated L-BFGS regularized inversion method replaces the traditional finite difference forward solver with an EFHPN forward model. This not only achieves millisecond-level acceleration in a single forward calculation but also directly calculates the gradient of the objective function using the automatic differentiation technique of deep learning frameworks, avoiding the explicit construction and storage of the Jacobian matrix. The inversion method specifically includes the following sub-steps:

[0143] Data Acquisition: Collect MT observation data for the target area;

[0144] Objective function construction: The objective function is a weighted combination of the data fitting term and the Tikhonov model regularization term. The data fitting term measures the difference between the observed data and the response data predicted by the neural operator forward model, while the model regularization term constrains the smoothness of the inversion results.

[0145] Iteration parameter settings: Set the maximum number of iterations and the iteration convergence criteria. For example, the maximum number of iterations can be set to 150, and a convergence criterion based on the gradient norm and the change in the objective function can be established.

[0146] Model update strategy: In each iteration, the L-BFGS algorithm updates the approximation of the inverse Hessian matrix through a double-loop recursive mechanism, and combines it with a line search that satisfies the strong Wolfe condition to determine the optimal step size, thereby completing the update of the model parameters until the convergence condition is met and the resistivity model of the inversion region is output.

[0147] Comparison of different forward proxy models:

[0148] The forward model constructed by the technical solution of this invention adopts a branch-backbone fusion structure, and introduces a hybrid periodic neural representation in the backbone network. To illustrate the technical effect of the model of this invention, this embodiment compares the key technical features of four MT forward surrogate models. All models use the same FNO six-layer branch network, and the differences are mainly reflected in the backbone network structure:

[0149] Table 1. Comparison of key technical features of different MT forward proxy models

[0150]

[0151] The backbone network of EFNO uses the conventional fully connected network FCNN and the GELU activation function, which has relatively limited ability to express rapidly changing or periodic oscillation features. The backbone network of EFPN uses a single periodic neural representation network, which can capture periodic oscillation features, but when using a periodic activation function alone, it may suffer from insufficient stability in some complex scenarios. The backbone network of EFKAN uses KAN, a hybrid activation of SiLU and B-spline functions, which has strong nonlinear fitting ability, but B-splines usually require grid discretization, resulting in high computational complexity. EFHPN is the preferred model proposed in this invention. Its backbone network uses a hybrid periodic activation mechanism of SiLU and SIREN, which does not require additional grid discretization, has moderate computational complexity, and has continuous coordinate prediction capabilities.

[0152] The HPN design in the EFHPN model is as follows: the SiLU branch is responsible for fitting the smooth decay trend of the electromagnetic response and the low-frequency background field, while the SIREN branch utilizes the inherent periodicity of the sine function to capture periodic oscillations and phase change characteristics. The outputs of the two branches are then element-wise summed and fused, and subsequently normalized using LayerNorm and activated by PReLU. This design improves the network's ability to represent complex electromagnetic responses while maintaining computational efficiency, and simultaneously preserves continuous prediction capability and training stability.

[0153] Dataset construction and network training strategies:

[0154] Deep learning model training heavily relies on high-quality "resistivity model-electromagnetic response" training samples. To endow forward models with good predictive power, the training set needs to cover multiple complex resistivity models that may exist in subsurface media. This invention employs a Gaussian random field generation strategy based on spectral methods, whose energy spectrum... Defined as:

[0155] ;

[0156] Where k represents the wavenumber and β is the scaling factor. By controlling the scaling factor, random fields with different spatial frequency characteristics can be generated. In this embodiment, random fields with five scales of β=[3,4,5,6,7] are used to construct a multi-scale resistivity model and generate a dataset containing 10,000 training samples and 1,000 test samples.

[0157] This invention further designs a forward computational domain grid partitioning strategy for generating training datasets. For example... Figure 3 As shown in (a), the complete two-dimensional computational domain spans 400 km horizontally and 300 km vertically, including a 100 km air layer and a 200 km underground extension region. Figure 3 As shown in (b), the target observation area uses a 128×64 structured grid. The core area has a uniform grid in the horizontal direction and a non-uniform grid in the depth direction that gradually increases downwards. The thickness of the first layer is set to 5m, and there are a total of 64 layers, increasing sequentially. The outer extension area and air layer use a 10-layer geometrically stretched grid, and the grid of the entire computational domain is 148×84. 64 observation points are uniformly distributed on the surface, with an observation frequency of 10... -3 Hz to 10 3 41 frequency points were selected at logarithmic intervals between Hz. An air layer was assigned 10... 9 The ultra-high resistivity value (Ω.m) is achieved by smoothly expanding outwards by five transition grids to the left, right, and downwards from the periphery of the core heterogeneous region. This, along with the surrounding stretched edge-expanding layer, is uniformly assigned a value of 10. 2 The modeling strategy employs a uniform background resistivity of Ω·m. This approach is beneficial for improving the numerical stability of the finite difference solution process.

[0158] Magnetotelluric observation data exhibits a wide frequency band distribution and cross-scale physical characteristics. Conventional absolute mean square errors may be dominated by large-amplitude responses, causing the network to ignore low-resistivity anomalies and phase details. To balance different response components, this invention employs a relative L2 loss function, which can be expressed as:

[0159] ;

[0160] In the formula, N is the number of training samples; and These are the apparent resistivity and phase reference solutions obtained by the finite difference method, respectively. and These are the apparent resistivity and phase predicted by the neural operator forward model, respectively.

[0161] The AdamW optimizer was used for network parameter optimization, with an initial learning rate set to 10. -3 The weight decay factor is 10. -4The learning rate is scheduled using a StepLR dynamic annealing strategy, decaying to 0.5 times its original value every 50 epochs. For example... Figure 4 As shown, EFHPN in the training set ( Figure 4 a) and test set ( Figure 4 b) in the above examples all show a stable convergence trend and the lowest global loss. Figure 4 (c) intuitively reflects the smooth decay process of the learning rate, effectively eliminating numerical oscillations in the later stages of training.

[0162] Theoretical model trial calculation:

[0163] This embodiment selects four network architectures—EFNO, EFPN, EFKAN, and EFHPN—for comparative calculations, and uses the numerical solution of the Finite Difference Method (FDM) as a benchmark to evaluate the predictive performance of the preferred model (EFHPN model) of this invention. To ensure consistency in the comparison, each deep learning model is trained on the same training set and uses consistent major hyperparameter configurations.

[0164] Figure 5 The prediction results of four neural operator forward models on the test samples are shown. Figure 5 (a) is the resistivity model of the test sample. Figure 5 The b1-b4 series of diagrams are the FDM reference solutions. Figure 5 The c1-c4 series represents the prediction results for EFNO, the d1-d4 series for EFPN, the e1-e4 series for EFKAN, and the f1-f4 series for EFHPN. Labels 1, 2, 3, and 4 correspond to four types of electromagnetic responses: TE apparent resistivity, TE phase, TM apparent resistivity, and TM phase, respectively. From the overall response morphology, all four models are able to capture the main spatial response characteristics.

[0165] To further compare the performance of different architectures, this embodiment uses relative error for apparent resistivity and absolute error for phase. The error calculation formula is as follows:

[0166] ;

[0167] in, and These represent the baseline solution of apparent resistivity calculated using the finite difference method and the apparent resistivity predicted by the neural operator forward model, respectively. and These are the phase reference solution calculated using finite difference and the phase predicted by the neural operator forward model, respectively. These are apparent resistivity error and phase error, respectively.

[0168] Figure 6The error distributions of four neural operator forward models on the test samples are shown, including the EFNO method. Figure 6 The a1-a4 values ​​exhibit significant bias in apparent resistivity prediction, with local relative errors exceeding 5% and phase errors of approximately 1-1.5 degrees. The EFPN method ( Figure 6 The prediction accuracy (b1-b4) is significantly improved by the introduction of a periodic sinusoidal activation function, reducing the relative error of apparent resistivity to within 2.5%. The EFKAN method ( Figure 6 The accuracy of (c1-c4) has also improved. In comparison, EFHPN ( Figure 6 The model (d1-d4) exhibits the best prediction performance. The results show that, compared to models using conventional activation functions, models employing periodic or mixed periodic representations can reduce errors in some regions. The error distribution of EFHPN is more uniform, indicating its better characterization ability for MT responses under complex resistivity structures.

[0169] One of the advantages of neural operators over conventional deep learning methods is their ability to predict arbitrary positions and frequencies within the training scope, rather than being limited to a fixed output grid during training. This capability stems from the backbone network's encoding mechanism for continuous coordinates, enabling the network to learn a continuous mapping relationship between coordinates and responses.

[0170] Figure 7 The results show the forward response predictions of four neural operator networks at two different measurement points (y = -15 km and y = 15 km). These locations are not the 64 discrete measurement points used during training, but arbitrary locations obtained through continuous spatial interpolation. The prediction results demonstrate that all four neural operator networks can accurately predict the electromagnetic response at any location within the measurement point range. EFHPN maintains excellent prediction accuracy at all measurement points, accurately capturing the spatial variation characteristics of the electromagnetic response and achieving a high degree of agreement with the finite difference FDM benchmark solution.

[0171] Figure 8This paper presents the forward response prediction results of four neural network operators along a survey profile at two specific frequency points (500Hz and 0.02Hz). These frequency points are not the 41 discrete frequency points used during training, but arbitrary frequencies obtained through continuous interpolation. The prediction results show that all neural operator methods can provide reasonable electromagnetic response predictions at arbitrary frequencies, closely matching the finite difference FDM benchmark solution. EFHPN exhibits superior prediction accuracy across all frequency bands, particularly in the high-frequency (500Hz) and low-frequency (0.02Hz) bands, where its prediction curve almost perfectly overlaps with the FDM solution. In contrast, other methods show slight deviations in certain frequency bands; for example, EFNO and EFKAN show subtle differences in apparent resistivity predictions at high frequencies. This validates the advantage of EFHPN's periodic activation function in handling broadband electromagnetic responses.

[0172] like Figure 9 As shown, this invention compares the acceleration of L-BFGS inversion of subsurface resistivity profiles when using four neural operators—EFNO, EFPN, EFKAN, and EFHPN—as forward modeling engines. Real model ( Figure 9 (a) It possesses the following geological features: a shallow high-resistivity layer, a middle transition zone, and a deep low-resistivity basement, with complex undulating geological interfaces, including the downward extension of the high-resistivity layer on the right and the obvious uplift of the low-resistivity basement in the middle (0-10 km). This model poses a significant challenge to the vertical resolution and morphological characterization capabilities of the inversion algorithm. Inversion results show that the inversion effects driven by different neural operators vary: the EFNO-accelerated L-BFGS inversion results blur the boundary between the deep low-resistivity basement and the overlying strata, failing to distinguish the undulating morphology. Figure 9 (b); The EFKAN-accelerated L-BFGS inversion results barely recovered the main layered outline, but distorted the true shape of the deep low-resistivity anomaly, flattened the originally undulating base, and had defects in detail depiction. Figure 9 c); EFPN accelerates L-BFGS inversion results, with better boundary resolution than the other two methods mentioned above, but still has shortcomings in local structural details (such as the connectivity and boundaries of the low-resistivity core region). Figure 9 (d); EFHPN accelerated L-BFGS inversion results showed the best underground model reconstruction effect, with the inverted profile highly consistent with the real model. It not only accurately reconstructed the spatial distribution of resistivity, but also captured the spatial undulation characteristics of the interface. In particular, the upward bulge of the low-resistivity base in the middle and the boundary of the high-resistivity anomaly on the right were clearly and sharply depicted, almost restoring all the structural details of the real model. Figure 9 (e).

[0173] The above comparison results demonstrate that the fundamental difference in inversion imaging resolution stems from the accuracy of different neural operators in forward mapping the electromagnetic field. EFHPN, with its hybrid periodic sinusoidal activation function design, can more sensitively capture electromagnetic field distortions caused by complex heterogeneous geological bodies, thus providing the L-BFGS optimization engine with high-confidence forward predictions and accurate physical gradient fields, guiding the algorithm to converge to the globally optimal solution that best approximates the true structure. This result confirms the decisive role of high-precision forward operators in the quality of nonlinear inversion.

[0174] In some embodiments, the present invention also provides a forward modeling system based on the above-described forward modeling method, including a model building module, a data acquisition module, and a forward modeling module that are connected in sequence or to each other.

[0175] The model building module is used to build a neural operator forward model that includes an FNO branch network, an HPN backbone network, and an operator fusion module.

[0176] The data acquisition module is used to acquire the underground two-dimensional resistivity model ρ(y,z) and the horizontal coordinate y and frequency f of the observation point;

[0177] The forward modeling module is used to input the subsurface two-dimensional resistivity model ρ(y,z) into the branch network of the forward model, and to input the horizontal coordinates y and frequency f of the observation points into the backbone network of the forward model. It then uses the operator fusion module to achieve a continuous mapping from the subsurface resistivity distribution to the surface electromagnetic response. The electromagnetic response includes the apparent resistivity in TE mode, the phase in TE mode, the apparent resistivity in TM mode, and the phase in TM mode.

[0178] In some embodiments, the present invention also provides an inversion system based on the above-described forward modeling method, comprising an inversion model building module, an observation module, and an inversion module that are connected sequentially or interconnected.

[0179] The inversion model building module is used to use the pre-trained neural operator forward model as a forward modeling engine to replace the numerical forward modeling solver.

[0180] The observation module is used to collect observation data of the inversion area through detection equipment;

[0181] The inversion module is used to perform L-BFGS iterative inversion based on observation data, the initial model, and the neural operator forward model to obtain the complex underground electrical structure of the inversion area;

[0182] In some embodiments, the system also includes a gradient calculation module for calculating the gradient of the objective function with respect to the model parameters using an automatic differentiation mechanism.

[0183] It should be understood that the specific implementation process of each module is described in the above method. This invention will not repeat the details here. The above division of functional modules is only for illustrative purposes. In some embodiments, some functional modules can be combined and some functional modules can be separated. Each functional module can be implemented in software, hardware, or a combination of software and hardware. The software and hardware devices include, but are not limited to, general-purpose computer equipment, programmable gate arrays, digital signal processors, microprocessors and their corresponding programming or burning software.

[0184] The present invention also provides a computer terminal, including one or more processors and a memory storing a computer program; wherein the processor calls the computer program to implement the steps of the above-mentioned neural operator-based magnetotelluric forward modeling method, or to implement the steps of the above-mentioned neural operator-based magnetotelluric inversion method.

[0185] Please refer to the explanation of the method described above for the specific implementation process of each step.

[0186] It should be understood that, in the embodiments of the present invention, the processor may be a Central Processing Unit (CPU), or it may be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor. The memory may include read-only memory and random access memory, and provides instructions and data to the processor. A portion of the memory may also include non-volatile random access memory. For example, the memory may also store device type information.

[0187] The present invention also provides a computer-readable storage medium storing a computer program, which, when called by a processor, implements the steps of the above-described magnetotelluric forward modeling method based on neural operators, or implements the steps of the above-described magnetotelluric inversion method based on neural operators.

[0188] Please refer to the explanation of the method described above for the specific implementation process of each step.

[0189] The computer-readable storage medium can be an internal storage unit of the hardware and software device in any of the foregoing embodiments (e.g., the hard drive or memory of the controller), or it can be an external storage device of the controller, such as a plug-in hard drive, smart memory card (SMC), secure digital card (SD card), flash memory card, etc., equipped on the controller. Furthermore, the storage medium can also simultaneously include both the controller's internal storage unit and external storage device.

[0190] Based on the above understanding, the core part of the technical solution of this invention that contributes to the prior art, or all or part of the content of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and contains several instructions to cause a computer device (such as a personal computer, server, or network device) to execute all or part of the steps of the methods described in the various embodiments of this invention. Available storage media include, but are not limited to: USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, optical disks, and other media capable of storing program code.

[0191] Those skilled in the art will understand that embodiments of this application can be provided in the form of a method, system, or computer program product. Therefore, this application can be implemented entirely in hardware, entirely in software, or a combination of hardware and software. Furthermore, this application can also be implemented as a computer program product containing computer-usable program code on a computer-readable storage medium (such as a disk storage device, CD-ROM, optical storage, etc.). Embodiments of this application are described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products. It should be understood that the function of each step in the flowchart and / or each block in the block diagram can be implemented by computer program instructions. These computer program instructions can be executed by a processor of a general-purpose computer, special-purpose computer, or other programmable data processing device to generate means for implementing the functions specified in the flowchart and / or block diagrams. These instructions can also be stored in a computer-readable storage medium to cause a computer or other programmable device to operate in a particular manner, thereby producing an article of manufacture containing instruction means to implement the functions specified in the flowchart and / or block diagram. Furthermore, the instructions can be loaded onto a computer or other programmable device to form a computer-implemented processing flow by performing a series of operational steps, so that the instructions executed on the computer or other programmable device can implement the functional steps specified in the flowchart and / or block diagram.

[0192] It should be emphasized that the examples described in this invention are illustrative rather than limiting. Therefore, this invention is not limited to the examples described in the specific embodiments. Any other embodiments derived by those skilled in the art based on the technical solutions of this invention, without departing from the spirit and scope of this invention, whether modifications or substitutions, are also within the protection scope of this invention.

Claims

1. A magnetotelluric forward modeling method based on neural operators, characterized in that: Includes the following steps: Step 1: Obtain the underground two-dimensional resistivity model ρ(y,z) and the horizontal coordinate y and frequency f of the observation point, where z represents the depth coordinate; Step 2: Input the data obtained in Step 1 into the neural operator forward model. The neural operator forward model adopts the branch-backbone dual network architecture of the DeepONet deep operator network, which includes a backbone network, branch networks, and an operator fusion module. The Fourier neural operator FNO is introduced into the branch network. The branch network is used to extract the global medium feature vector b from the two-dimensional resistivity model ρ(y,z) of the subsurface. k k is the feature channel index; the backbone network introduces a hybrid periodic neural representation dual-branch fusion mechanism, and the backbone network is used to generate continuous coordinate basis functions t based on the continuous observation coordinates (y,f). k ; The operator fusion module will integrate the continuous coordinate basis function t k With global medium eigenvector b k By fusing tensor inner products, the electromagnetic response at the corresponding observation location and frequency is output, thus obtaining the forward modeling result; The backbone network comprises a normalization layer, multiple hidden fully connected layers, and an output layer connected in sequence. The normalization layer normalizes the horizontal coordinate y and frequency f of the observation points. The multiple hidden fully connected layers perform high-dimensional nonlinear feature mapping on the normalized continuous observation coordinates (y, f). The hidden fully connected layers employ a hybrid periodic neural representation dual-branch fusion mechanism, where the two branches are a smooth branch and a sinusoidal periodic branch. The sinusoidal periodic branch uses a periodic activation function, and the smooth branch uses a nonlinear activation function. The output layer generates the continuous coordinate basis function t. k； The electromagnetic response includes TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase.

2. The method according to claim 1, characterized in that: The branch network has a sequentially connected up-dimensional layer P, several cascaded Fourier neural operator layers, and a down-dimensional layer Q; Each Fourier neural operator layer employs the Fourier neural operator FNO, comprising two parallel branches. The first branch follows the sequence: "Fourier transform F → frequency domain spectral weighting modulation R → inverse Fourier transform F". -1 The frequency domain convolution paradigm is used; the second branch uses a local linear transformation W as a bypass connection; finally, the outputs of the two branches are summed and fused to obtain the latent feature tensor. The latent feature tensor obtained from the current layer is used as the input latent feature tensor of the next Fourier neural operator layer, and the high-dimensional latent feature tensor output from the last Fourier neural operator layer is then input to the dimensionality reduction layer Q.

3. The method according to claim 1, characterized in that: The smooth branch uses the SiLU activation function, and the sinusoidal periodic branch uses the sinusoidal activation function. The outputs of the two branches are summed and fused element-wise. The corresponding hybrid periodic neural representation is as follows: ； in, This is a branch of the SiLU activation function, used to characterize the smooth decay trend and background changes in the electromagnetic response; W1 and W2 are the sinusoidal activation function branches used to characterize the periodic oscillations and phase changes in the electromagnetic response; W1 and W2 are learnable weight matrices, and ω0 is a hyperparameter controlling the frequency scale of the sinusoidal function. denoted as the hybrid activation function in the hidden fully connected layer of the backbone network; x represents the input of the hidden fully connected layer.

4. A magnetotelluric inversion method based on neural operators, characterized in that: Includes the following steps: Construct an inversion framework based on a neural operator forward model to accelerate the inversion process: use the pre-trained neural operator forward model as described in any one of claims 1-3 as the forward modeling engine; Define the objective function and the initial model parameters m0; Acquisition and preprocessing of observation data: Observation data d of the inversion area were acquired using MT detection equipment. obs ; Iterative inversion: based on initial model parameters m0, a neural operator forward model as the forward modeling engine, and observation data d. obs Perform iterative inversion and update the model parameters m; Here, the model parameter m is the target parameter model of the inversion area, which represents the complex underground electrical structure; the model parameter m obtained at the end of the iteration is used as the inversion result of the underground electrical structure of the inversion area.

5. The method according to claim 4, characterized in that: The inversion process employs L-BFGS regularized inversion, and the objective function is constructed using the Tikhonov regularization framework. The L-BFGS regularized inversion process is as follows: Based on the initial model parameters m0, the neural operator forward model is invoked to calculate the model parameters m corresponding to the k-th iteration. k The corresponding predicted response; Based on predicted response and observation data d obs Construct the objective function And use automatic differentiation to calculate the objective function. Relative to m k gradient k is the iteration round marker; Gradient based on the objective function The L-BFGS algorithm is used to update the search direction p. k The step size α is determined by line search. k This updates the model parameters corresponding to the (k+1)th iteration. Enter the iterative operation; The iteration terminates when the gradient norm is less than a preset threshold or the change in the objective function is less than a preset threshold. The model parameters m at the end of the iteration inv Output as the target parameters of the inversion region.

6. A forward modeling system based on the method of any one of claims 1-3, characterized in that: include: The model building module is used to build a neural operator forward model that includes a branch network, a backbone network, and an operator fusion module. The data acquisition module is used to acquire the underground two-dimensional resistivity model ρ(y,z) and the horizontal coordinate y and frequency f of the observation point; The forward modeling module is used to input the subsurface two-dimensional resistivity model ρ(y,z) into the branch network of the neural operator forward model, and input the horizontal coordinate y and frequency f of the observation point into the backbone network of the neural operator forward model. The operator fusion module realizes the continuous mapping from the subsurface resistivity distribution to the surface electromagnetic response. The electromagnetic response includes TE mode apparent resistivity, TE mode phase, TM mode apparent resistivity, and TM mode phase.

7. An inversion system based on the method of any one of claims 4-5, characterized in that: include: The inversion model building module is used to use the pre-trained neural operator forward model as the forward modeling engine; The observation module is used to collect observation data of the inversion area through the MT detection equipment. obs ; The inversion module is used to perform forward modeling based on initial model parameters m0, a neural operator forward model serving as the forward modeling engine, and observation data d. obs The model parameters m are updated through iterative inversion to finally obtain the complex underground electrical structure of the inverted area.

8. A computer terminal, characterized in that: include: One or more processors; A memory that stores computer programs; The processor invokes a computer program to achieve the following: The steps of a magnetotelluric forward modeling method based on neural operators as described in any one of claims 1-3 or the steps of a magnetotelluric inversion method based on neural operators as described in any one of claims 4-5.

9. A computer-readable storage medium, characterized in that: The computer program is stored and is invoked by the processor to implement: The steps of a magnetotelluric forward modeling method based on neural operators as described in any one of claims 1-3 or the steps of a magnetotelluric inversion method based on neural operators as described in any one of claims 4-5.

Images