A method for neural operator network training and modeling for wave equation solving

Abstract

This disclosure provides a self-supervised training method for neural operator networks used to solve wave equations. The method includes the following steps: combining wave velocity field, source location, and frequencies obtained through dynamic random sampling into an input sample; selecting a set of input samples from the dataset to form a training batch; preprocessing the training batch data to a suitable format and inputting it into the neural operator network; calculating the forward propagation of the neural network and outputting a predicted wave field; using the predicted wave field as the initial field for Black-Scholes (BS) iterations, performing multiple BS iterations based on a dynamic iteration step number k, and outputting multiple residual fields; taking the L2 norm of the k residual fields and summing them, outputting the result, and determining whether the loss has reached the convergence condition; if the convergence condition is reached, exiting the loop and outputting the trained neural operator network. This avoids the cost of preparing labeled data in existing supervised training methods and overcomes the convergence problem of existing BS methods.

Description

Technical Field

[0001] This disclosure relates to the field of computer science, and more specifically, to a method for training and modeling neural operator networks for solving wave equations. Background Technology

[0002] Wave equations are a collective term for equations that describe wave phenomena, including the equations for sound waves, elastic waves, and electromagnetic waves.

[0003] A wave field is a physical field formed by the distribution of wave physical quantities (such as pressure) in spacetime. The goal of solving the wave equation is to obtain the wave field distribution.

[0004] A neural operator (or neural operator network) is a special type of neural network used to process mappings from physical fields to physical fields.

[0005] Wave phenomena described by wave equations are prevalent in scientific and engineering fields, and wave equation solvers have wide applications. Wave equations include those for sound waves, elastic waves, and electromagnetic waves, and are used in fields such as medical ultrasound, geological exploration, and electrical equipment design. Solving wave equations is of great significance in engineering. In seismic exploration, seismic waves are generated from the Earth's surface / sea surface. These waves propagate underground, undergo reflection and refraction, and return to the Earth's surface / sea surface, where they are received by a receiver. Using the received reflected and refracted wave signals, the subsurface structure can be imaged using full waveform inversion technology, thereby determining oil and gas reserves.

[0006] Wave equation solvers generally employ frequency domain and time domain algorithms. Frequency domain algorithms can avoid errors caused by time discretization and iterative calculations with time progression, making them suitable for analyzing specific frequency components and thus widely used. Among frequency domain algorithms, the Born Series method offers high numerical accuracy, low memory requirements, and ease of adaptation to different computing hardware architectures. The open-source acoustic wave simulator j-Wave provides a Born Series solver implementation based on Google JAX. Among existing technologies, the Born Series (BS) method is a high-precision method for solving wave equations. Its idea is to introduce a scattering potential variable to divide the non-homogeneous parameter field into a homogeneous part and a scattering part; the scattering part is considered as the source term for solving the homogeneous parameter field, and the true scattering solution is gradually approximated through iteration. It can solve various forms of wave equations. Its advantages are: high numerical accuracy, employing spectral methods in all directions, theoretically achieving arbitrary accuracy; and low spatial complexity. It has low memory requirements and low dependence on complex operators, making it easy to adapt to different computing hardware architectures (the main operators are FFT and vector multiplication). Its limitations are: it is applicable to a limited range of working conditions, and the series diverges when scattering is strong, making it impossible to compute. Therefore, its practical application is limited to weak scattering scenarios.

[0007] Besides traditional black-and-white (BS) methods, some studies have proposed combining AI with BS. The general idea behind these methods is to use neural networks to establish a mapping between the input parameters and output solutions of the equation, with the input parameter field and initial unsolved variables, and the output unsolved variable field. Supervised training is performed using a pre-constructed "parameter field – unsolved variable field" dataset. The advantage of this type of method is its fast inference speed. Its limitations are: reliance on training data, requiring a trade-off between the richness of the data distribution and the cost of data preparation; and existing research only solves standard Helmholtz-form univariate acoustic wave equations.

[0008] Existing Black-Scholes (BS) methods have limited applicability to certain conditions; when scattering is strong, the series diverges and cannot be calculated, thus limiting their practical application to weak scattering scenarios. Current AI methods combined with BS rely on training data, requiring a trade-off between the richness of the data distribution and the cost of data preparation; furthermore, current research only solves the standard Helmholtz form of univariate acoustic wave equations. Summary of the Invention

[0009] The purpose of this disclosure is to provide a self-supervised training method for neural operator networks for solving wave equations. It aims to construct a self-supervised physical loss function using the intermediate iteration process of the BS solver and dynamically adjust the parameter sampling range and the number of iterations, so as to achieve generalization model training at a lower overall cost than the prior art. This solves the problems of series non-convergence and the trade-off between data richness and data preparation cost in the prior art.

[0010] In general, a self-supervised training method for neural operator networks for solving wave equations is provided, comprising the following steps: Step 1: Define and initialize the neural operator network for processing the probed wave velocity field of matter or object; construct the wave velocity field dataset and determine the source location and frequency to be calculated for each wave velocity field; determine the dynamic iteration step k of BS in physics loss, and the computational parameters required for BS. Step 2: Combine a wave velocity field, a seismic source location, and a frequency obtained through dynamic random sampling into an input sample, and take a set of input samples from the dataset to form a training batch. Step 3: Preprocess the training batch data to make it into a suitable format and input it into the neural operator network; Step 4: Calculate the forward propagation of the neural network and output the predicted wave field; Step 5: Using the predicted wave field as the initial field for BS iteration, perform multiple BS iterations based on the dynamic iteration step number k to output multiple residual fields. Step 6: Take the L2 norm of the k residual fields and sum them to output the physics loss value; Step 7: Determine whether the physics loss has reached the convergence condition. If it has, exit the loop and output the trained neural operator network. If it has not reached the convergence condition, proceed to step 8. Step 8: Call the automatic differentiation function provided by the programming framework to calculate backpropagation and obtain the gradient of the neural operator network parameters; Step 9: Call the optimizer to update the parameters of the neural operator network, and return to step 2 to repeat the process.

[0011] The specific method for obtaining the frequency through dynamic random sampling is as follows: the training is divided into several stages, and a different frequency sampling range is set for each stage. When the training in a certain stage converges, the next stage is entered.

[0012] The specific method for multiple BS iterations based on dynamic iteration step k is as follows: the training is divided into several stages, and the BS iteration step k of the physical loss function is dynamically adjusted with each training round. Different k values ​​are set for each stage. When the training in a certain stage converges, the next stage of training begins.

[0013] The neural operator network is structured as follows: it adopts an axially adaptive Fourier neural operator architecture, specifically including: The input data first passes through an input smoothing projection layer, then through a feature block and position encoding module, before entering the input smoothing projection layer. The core backbone network is composed of layer-axial-AFNO modules; The axial-AFNO module introduces a shared projection substrate, shared across different spatial axes in the input data. Simultaneously, a specific gating mechanism is introduced for each axis. Features, after frequency domain modulation, are restored via inverse Fourier transform. The specific implementation of the shared projection substrate in the axial-AFNO module is as follows: For the multidimensional features of the input, the network constructs a globally shared, learnable projective weight matrix. ; Specifically, after the input features undergo a one-dimensional real-valued Fast Fourier Transform along a specific spatial axis, they are divided along the channel dimension. Each group contains 1 independent group. There are 1 channel; for the extracted low-frequency component, let its feature within a given group be . Where L is the number of cutoff frequency modes, the shared basis is linearly reduced and increased in dimensionality through low-rank projection; the features are separated into real parts. and the virtual part Subsequently, shared matrices were used. , Given a set low-rank size, and k as the current layer number, perform spatial orientation-independent projection calculations: After subsequent gated modulation, a shared matrix is ​​then used. Map the rank dimension back to the feature dimension.

[0014] The specific implementation of the aforementioned gating mechanism is as follows: A specific gating mechanism is introduced for each axis, specifically by assigning specific, learnable gating parameters and biases to each independent spatial axis in the frequency domain. Furthermore, this gating mechanism introduces a frequency band partitioning strategy. (The truncated portion...) The low-frequency mode is divided into several frequency bands according to frequency level, for the first... The b-th frequency band of each spatial axis is defined with a specific gate vector as follows: After passing through the shared base After projection, the specific gate vector is used to perform element-wise multiplication modulation with the feature: The modulated features are processed by a nonlinear activation function and the above projection and gating process is repeated before entering the next layer; Finally, the frequency domain features processed independently on each axis are restored to the spatial domain by inverse discrete Fourier transform. When multi-axis feature fusion is used, learnable fusion weight parameters are initialized for each spatial axis. After normalization by the Softmax function, the output wavefields of each axis are weighted and summed to obtain the final wavefield prediction result that balances computational efficiency and anisotropic features.

[0015] The smooth projection layer contains stacked convolutional kernels with sizes of [sizes not specified]. 3D convolution.

[0016] In another general aspect, an end-to-end full waveform modeling and updating method is provided, comprising four processing stages: Phase 1: Physics-driven input construction, starting with the velocity model of the input. Random sampling and enhancement are performed, followed by combining the angular frequencies of the random samples. and the location of the epicenter Perform dimensionless processing to generate a normalized velocity model. and normalized source mask ; Phase Two: Parameter-efficient neural prediction. The normalized input first passes through a pre-computation module, which analytically calculates a dimensionless background wave field. To address source singularities, the input data is then fed into a frequency band partitioning expert system. This system comprises multiple parallel expert sub-models, each outputting a predicted value for the scattered wavefield. Finally, the background wavefield is superimposed to obtain the full wavefield prediction. ; Phase 3: CBS-guided physical supervision, constructing an unlabeled physical loss function, and converting the predicted wave field output by the network into a physical loss function. As initial values, a convergent Born series is iterated with a fixed number of steps, and the CBS operator is used. Using the fixed-point property of the wave field, calculate the consistency error between the predicted wave field and the iterated wave field; Phase Four: Velocity inversion based on reflection sensing. The parameters of the neural operator network trained using a self-supervised training method for solving wave equations are frozen and used as a differentiable forward modeling engine, introducing the up-wave separation operator. When calculating the backpropagation gradient, the energy of the uplink reflected wave is explicitly separated and weighted, thereby enabling the recovery of the fine underground structure from strong direct wave interference when updating the velocity model.

[0017] The technical effects to be achieved by the embodiments of the present invention are as follows: (1) The output of the finite steps in the Born Series iteration process is used to construct the physical loss function for self-supervised training of the neural network, thereby avoiding the cost of preparing labeled data in the existing supervised training method and overcoming the non-convergence problem of the existing BS method.

[0018] (2) During the implementation of the training method, parameters such as earthquake source location and frequency are sampled online in a dynamic and random manner to reduce the difficulty of training the neural operator network.

[0019] (3) During the above training process, the number of iterations of the physical loss function is dynamically adjusted with each training round to reduce the training difficulty of the neural operator network and accelerate the training process. Attached Figure Description

[0020] The above and other objects and features of this disclosure will become clearer from the following description taken in conjunction with the accompanying drawings.

[0021] Figure 1 This is a schematic diagram illustrating an example of a Born Series method flowchart in the prior art; Figure 2 This is a schematic diagram illustrating an example of the flowchart for combining AI with Born Series in existing technologies; Figure 3 This is a schematic diagram illustrating the training and inference process of a neural operator network based on Born series (BS) physical loss according to this disclosure; Figure 4This is a schematic diagram illustrating an end-to-end full waveform modeling and inversion framework architecture according to the present disclosure; Figure 5 This is a schematic diagram illustrating the internal architecture of a neural operator network according to this disclosure; Figure 6 This is a schematic diagram illustrating the architecture of a self-supervised training method for solving wave equations using neural operator networks according to Embodiment 1 of this disclosure; Figure 7 This is a schematic diagram showing wavefields at different source locations and frequencies according to Embodiment 2 of this disclosure; Figure 8 This is a schematic diagram illustrating the architecture of a self-supervised training method for solving wave equations using a neural operator network according to Embodiment 2 of this disclosure; Figure 9 This is a schematic diagram illustrating the 2D frequency domain forward modeling accuracy verification effect under complex media according to Embodiment 3 of this disclosure; Figure 10 This is a schematic diagram illustrating the effect of full waveform inversion application of out-of-distribution data in Embodiment 4 according to the present disclosure; Figure 11 This is a schematic diagram illustrating the effect of scalability verification of a large-scale 3D forward modeling simulation according to Embodiment 5 of this disclosure. Detailed Implementation

[0022] The following detailed embodiments are provided to aid the reader in gaining a comprehensive understanding of the methods, apparatus, and / or systems described herein. However, various changes, modifications, and equivalents of the methods, apparatus, and / or systems described herein will become apparent upon understanding this disclosure. For example, the order of operations described herein is merely illustrative and is not limited to those orders set forth herein, but may be changed as will become clear upon understanding this disclosure, except for operations that must occur in a specific order. Furthermore, for clarity and conciseness, descriptions of features known in the art may be omitted.

[0023] The features described herein may be implemented in different forms and should not be construed as limited to the examples described herein. Rather, the examples described herein are provided only to illustrate some of the many feasible ways of implementing the methods, apparatus, and / or systems described herein, which will become clear upon understanding the disclosure of this application.

[0024] As used herein, the term “and / or” includes any one of the associated listed items and any combination of any two or more.

[0025] Although terms such as “first,” “second,” and “third” may be used herein to describe various components, assemblies, regions, layers, or parts, these components, assemblies, regions, layers, or parts should not be limited by these terms. Rather, these terms are used only to distinguish one component, assembly, region, layer, or part from another. Thus, without departing from the teaching of the examples described herein, the first component, first assembly, first region, first layer, or first part referred to as the first component, first assembly, first region, first layer, or first part may also be referred to as the second component, second assembly, second region, second layer, or second part.

[0026] In the specification, when an element (such as a layer, region, or substrate) is described as being "on" another element, "connected to," or "bonded to" another element, the element may be directly "on" another element, directly "connected to," or "bonded to" the other element, or one or more other elements may be present in between. Conversely, when an element is described as being "directly on" another element, "directly connected to," or "directly bonded to" another element, no other elements may be present in between.

[0027] The terminology used herein is for the purpose of describing various examples only and is not intended to limit disclosure. Unless the context clearly indicates otherwise, the singular form is intended to include the plural form as well. The terms “comprising,” “including,” and “having” indicate the presence of the described features, quantities, operations, components, elements, and / or combinations thereof, but do not preclude the presence or addition of one or more other features, quantities, operations, components, elements, and / or combinations thereof.

[0028] Unless otherwise defined, all terms used herein (including technical and scientific terms) shall have the same meaning as commonly understood by one of ordinary skill in the art to which this disclosure pertains upon understanding this disclosure. Unless expressly defined herein, terms (such as those defined in a general dictionary) shall be interpreted as having a meaning consistent with their meaning in the context of the relevant field and in this disclosure, and shall not be interpreted in an idealized or overly formalistic manner.

[0029] Furthermore, in the description of the examples, detailed descriptions of well-known related structures or functions will be omitted when it is believed that such detailed descriptions would lead to a vague interpretation of this disclosure.

[0030] This invention aims to enable the detection of matter or objects in the field of geophysics using specific solution methods. It is applicable to scenarios where AI4S solves wave equations, including full waveform inversion in seismic exploration and medical ultrasound imaging. It can be integrated as a software module into software systems that rely on wave equation solvers, such as those used in seismic exploration and brain ultrasound CT imaging.

[0031] The frequency domain expressions of wave equations solvable by the BS method can all be written in a unified form. in The variable field to be determined is one or more physical fields with spatial and temporal distributions, characterizing wave physical quantities such as sound pressure and electric field strength; A parametric field is one or more physical fields with spatial or temporal distributions, which generally characterize the properties of the wave propagation medium. and It is a function related to a specific equation. It is a spatial differential operator containing only spatial uniformity coefficients, and its value depends on the specific equation. Form and All of these are known. It is unknown; the input to the wave equation solver is... The output is .

[0032] The solution approach and process of the BS method are as follows. Parameter field It is spatially non-uniform, meaning there is some scattering. The parametric field is divided into homogeneous components. (Unchanging with time and space) and the scattering component, thus the equation can be written as in This is called the scattering potential. To control the stability of the iteration, a small imaginary quantity is introduced. An iterative series can be constructed for the above equation. To gradually approach the final solution : The right-hand side of an iterative expression can be abstracted as the input. Output The iterative module BS Block, the entire iterative solution process is as follows: Figure 1 As shown.

[0033] The reason existing Black-Scholes (BS) methods fail to converge is that the spectral norm of the iterative operator is greater than 1. Multiple combinations cause the overall spectral norm to increase exponentially, eventually exceeding the machine's precision and resulting in NaN. Therefore, there are two ways to avoid computational divergence: one is to reduce the spectral norm of the iterative operator to less than 1; the other is to reduce the number of iterations and exit the iteration before divergence occurs. This invention adopts the second approach, using only a small number of iterations (reduced from thousands to less than 20, with the specific number of steps being an adjustable parameter) to construct the physics loss, thereby avoiding the problem of computational failure due to series non-convergence.

[0034] The principle of the overall scheme of this invention is as follows: Figure 3As shown, the principle of constructing physics loss through a small number of iterations in the BS iteration process is as follows. If the iteration level... For a convergent series, then hour Before the series converges, This actually represents the residual of the equation. The Black-Scholes iteration process is a process of gradually reducing the residual of the equation. When the residual is 0, it represents... Satisfy the equation, that is It is the solution to the equation. (This invention only considers the case where the equation has a unique solution and is within a smooth neighborhood. Since the physical process described by the wave equation is always smooth and bounded, this condition applies to most cases.) The basis for constructing the physics loss in this invention is the equation residual. When the loss is 0, the equation residual is also 0. Therefore, optimizing this loss can enable the neural network to learn the solution to the equation.

[0035] The specific end-to-end full waveform modeling and inversion framework is as follows: Figure 4 As shown, unlike traditional numerical iterative solutions, this framework includes the following four core processing stages: Phase 1: Physics-Driven Input Construction. To decouple the physical scale and enhance generalization, the velocity model of the input is first constructed. Random sampling and enhancement are performed. Then, the angular frequency of the random sampling is combined. and the location of the epicenter Perform dimensionless processing to generate a normalized velocity model. and normalized source mask This step eliminates the explicit dependence of the wave equation on specific frequency values, enabling the neural operator network to learn wave field features within a unified normalized space.

[0036] Phase Two: Parameter-Efficient Neural Prediction. The normalized input first passes through a pre-computation module, which analytically calculates a dimensionless background wave field. To address source singularities, the input data is then fed into a frequency band partitioning expert system. To resolve spectral bias issues in broadband modeling, this system comprises multiple parallel expert sub-models (Expert 1 to Expert M), each configured to focus on a specific narrow frequency band (e.g., dividing 4-20Hz into multiple sub-bands). These expert sub-models output predicted values ​​of the scattered wavefield, which are then superimposed with the background wavefield to obtain a full-field prediction. .

[0037] Phase 3: CBS-guided Physical Supervision. This phase constructs an unlabeled physical loss function. Unlike traditional solvers that iterate until convergence, this method uses the predicted wavefield output by the network. As an initial value, expand a fixed number of steps (e.g.) The convergent Born series (CBS) iteration is performed (step by step). The CBS operator is used... The fixed-point property is used to calculate the consistency error between the predicted wavefield and the iterated wavefield. This process utilizes the constraints of the physical equations themselves to guide the updating of the neural operator network parameters, without requiring real wavefield labels.

[0038] Phase Four: Velocity Inversion Based on Reflection Sensing. In the inversion phase, the parameters of the trained neural operator network are frozen and used as a differentiable forward modeling engine. An up-wave separation operator is introduced. When calculating the backpropagation gradient, the energy of the uplink reflected wave is explicitly separated and weighted, thereby enabling the recovery of the fine underground structure from strong direct wave interference when updating the velocity model.

[0039] Internal architecture of neural operator networks: To further improve the computational efficiency and feature extraction capability of the model when processing large-scale data, the neural operator network in this embodiment adopts the following... Figure 5 The axial adaptive Fourier neural operator (axial-AFNO) architecture is shown. This is the key innovation that distinguishes this invention from general convolutional neural networks or conventional Transformers: Hierarchical expert structure: Each frequency band expert sub-model adopts a hierarchical structure similar to U-Net. Input data first passes through an input smoothing projection layer (containing stacked convolutional kernels with sizes of [missing information]). 3D convolutions are used to suppress mesh artifacts. This is followed by patching and positional encoding before entering the system... The core backbone network is composed of layer-axial-AFNO modules.

[0040] Frequency Domain Gating and Feature Hybridization: In the frequency domain, this method does not use a fully connected dense matrix, but instead introduces a shared-basis projection. These basis matrices are shared across different spatial axes to enforce physical consistency and reduce the number of parameters. Simultaneously, specific gates are introduced for each axis to capture the directional anisotropy of seismic wave propagation. Features are frequency-domain modulated and then reconstructed using inverse Fourier transform (iFFT). This design reduces computational complexity to [missing information]. The level enables the model to process large-scale 3D data with limited video memory.

[0041] The specific implementation method of the shared projection substrate is as follows: The specific implementation method of introducing a shared projection substrate into the axial-AFNO module is as follows: For multidimensional features of the input, the network does not perform operations on each independent spatial axis (such as two-dimensional axes). Axial or three-dimensional Instead of assigning independent weights to axes, a globally shared, learnable projective weight matrix is ​​constructed. .

[0042] Specifically, after the input features undergo a one-dimensional real-valued Fast Fourier Transform (rFFT) along a specific spatial axis, they are divided along the channel dimension. Each group contains [number] independent blocks, and each block contains [number] blocks. There are 10 channels. For the extracted low-frequency component, let its feature within a given group be _____. (Where $L$ represents the number of truncated frequency modes), the shared basis is linearly reduced and increased in dimensionality through low-rank projection. Taking the first-layer projection as an example, the features are separated into real parts. and the virtual part Subsequently, shared matrices were used. ( Perform spatial orientation-independent projection calculations for a given low-rank size: After subsequent gated modulation, a shared matrix is ​​then used. Mapping the rank dimension back to the feature dimension. This sharing mechanism forces the model to extract physically consistent wave features across different wave propagation directions, avoiding dense three-dimensional complex matrix multiplication and decoupling the parameters of the spatial dimension. The specific implementation of the gating mechanism is as follows: a specific gating mechanism is introduced for each axis. Specifically, although the projection base is globally shared, to capture the directional anisotropy of seismic wave propagation in different dimensions, the module assigns specific, learnable gate parameters and biases to each independent spatial axis in the frequency domain. Furthermore, this gating mechanism introduces a band index strategy. (The truncated text is incomplete and requires further context.) The low-frequency mode is divided into several frequency bands (Num Bands) according to frequency level. The $b$-th frequency band of the $n$ spatial axes is defined with a specific gate vector as follows: After passing through the shared base After projection, the specific gate vector is modulated by element-wise multiplication (Hadamard product) with the features: The modulated features are passed through a nonlinear activation function (such as GELU) and the above projection and gating process is repeated before entering the second layer.

[0043] Finally, the frequency domain features processed independently along each axis are restored to the spatial domain using the inverse discrete Fourier transform (iFFT). When multi-axis feature fusion is employed, learnable fusion weight parameters are initialized for each spatial axis. After normalization using the Softmax function, the output wavefields of each axis are weighted and summed to obtain a final wavefield prediction result that balances computational efficiency and anisotropic characteristics. Example 1

[0044] Figure 6 This is a schematic diagram illustrating a self-supervised training method for neural operator networks used to solve wave equations according to Embodiment 1 of this disclosure.

[0045] Taking the acoustic wave equation as an example, this invention is applied to train a neural network.

[0046] Step 1: Define and initialize the neural operator network; construct the wave velocity field dataset and determine which source locations and frequencies each wave velocity field needs to be calculated for; determine the number of iteration steps k of BS in the physics loss, as well as other computational parameters required by BS (if applicable).

[0047] Step 2: A wave velocity field, a seismic source location, and a frequency are combined to form an input sample. A set of input samples is taken from the dataset to form a training batch.

[0048] Step 3: Process the training batch data into a suitable format and input it into the neural operator network.

[0049] Step 4: Calculate the forward propagation of the neural network and output the predicted wave field.

[0050] Step 5: Use the predicted wave field as the initial field for the BS iteration, perform k BS iterations, and output k residual fields.

[0051] Step 6: Take the L2 norm of the k residual fields and sum them to output the loss value.

[0052] Step 7: Determine if the loss condition has been met. If it has, exit the loop and output the trained neural operator network.

[0053] Step 8: Utilize the automatic differentiation function provided by the programming framework to calculate backpropagation and obtain the gradients of the neural operator network parameters. (Note that this step requires both the neural operator network and the Black-Scholes solver used in the loss calculation to be implemented using a differentiable programming framework.) Step 9: Call the optimizer to update the parameters of the neural operator network. Return to step 2 and repeat. Example 2

[0054] Improvements were made to the input parameter selection strategy and the number of BS iteration steps (mainly steps 2 and 5) in Implementation Example 1.

[0055] Taking the equation of sound waves as an example, Figure 7 This diagram illustrates wavefields at different source locations and frequencies. The input parameter c includes the velocity field, source location, and frequency. Frequency determines the scale characteristics of the wavefield; the wider the frequency distribution, the wider the scale distribution of the wavefield, and the more difficult the corresponding learning task. Therefore, in the early stages of training, only a narrow frequency range can be sampled; as training progresses, the frequency range is gradually widened. This invention calls this parameter sampling method "dynamic random sampling." The specific strategy is to divide the training into several stages, each with a different frequency sampling range. When the training in a certain stage converges, the next stage begins. For example, if the final model is expected to have interpolation capability in the 5-10Hz frequency band, the training is divided into 5 stages. Phase 1: Randomly sample frequencies within the 5-6Hz range and train until convergence; Second stage: Randomly sample frequencies within the range of 5~7Hz and train until convergence; The third stage: Randomly sample frequencies within the range of 5~8Hz and train until convergence; Fourth stage: Randomly sample frequencies within the range of 5~9Hz and train until convergence; Fifth stage: Randomly sample frequencies within the range of 5~10Hz and train until convergence.

[0056] Another improvement of this invention is that, during training, the number of Black-Scholes iterations (k) of the physical loss function is dynamically adjusted with each training round. A smaller k results in a simpler loss function, easier optimization, and faster computation; a larger k amplifies the error in the predicted wave field, forcing the model to output a more accurate predicted wave field, but also makes the optimization problem more complex and computationally time-consuming. Therefore, a smaller k value can be used in the early stages of training, gradually increasing it as training progresses. Specifically, the training is divided into several stages, each with a different k value. When a stage converges, the next stage begins. For example, the training process can be divided into 5 stages. Phase 1: Set k = 5 and train until convergence; Second stage: Take k = 10 and train until convergence; Phase 3: Set k = 15 and train until convergence; Phase 4: Set k = 20 and train until convergence; Fifth stage: Take k = 30 and train until convergence.

[0057] Integrating the above improvements into the process of Embodiment 1 yields Embodiment 2, as follows: Figure 8 As shown. Example 3

[0058] 2D frequency domain forward modeling accuracy verification in complex media This embodiment verifies the accuracy and generalization ability of the neural operator network proposed in this disclosure in 2D frequency domain forward modeling through comparative experiments on public datasets and high-difficulty geological models.

[0059] Experimental setup: The OpenFWI benchmark dataset was selected for training, and the geologically complex Marmousi model was introduced as the out-of-distribution (OOD) test object. Comparison benchmarks included the traditional numerical solver (CBS), Physical Information Neural Network (PINN), Physically Augmented Neural Operator (PINO), and Data-Driven Operator Learning (DPOT).

[0060] Results analysis: Numerical precision: such as Figure 9 As shown, the method disclosed in this publication (CBNO) exhibits extremely high reconstruction accuracy on complex Marmousi models. Compared to the physically augmented PINO method, the root mean square error (RMSE) of this method is reduced by approximately 50.00% (from 0.0018 to 0.0009), and the relative error (RelErr) is reduced from 0.2707 to 0.1464.

[0061] Structural fidelity: Our method achieved 0.940 and 0.986 in structural similarity (SSIM) and Pearson correlation coefficient (Corr) respectively, which are significantly better than PINN (SSIM is only 0.395).

[0062] Frequency domain details: Wave field snapshots show that the wave field texture generated by this method is clear and accurately captures the complex interference fringes caused by faults and abrupt velocity changes, overcoming the "spectral deviation" problem of existing methods in broadband. Example 4

[0063] Full waveform inversion application of out-of-distribution (OOD) data This embodiment verifies the effect of using the neural operator network trained in this disclosure as a forward modeling engine in 2D full waveform inversion (FWI).

[0064] Inversion process: The parameters of the pre-trained neural operator network are frozen and embedded into the inversion loop. An up-wave separation operator is introduced during gradient calculation to extract the up-wave component carrying reflection information to construct the loss function, thereby improving deep imaging capabilities.

[0065] Results analysis: Generalization ability: As shown in Figure 10, purely data-driven methods (such as InverseNet and VeloGAN) suffer severe performance collapse when faced with the unseen Marmousi model, resulting in the loss of geological structure in the inversion results. In contrast, the method disclosed in this paper successfully reconstructs the complex layered structure and fault details of the Marmousi model.

[0066] Quantitative metrics: In the Marmousi test, the structural similarity (SSIM) of the inversion results of this method to the real velocity model reached 0.515, which is much higher than InverseNet's 0.310 and second only to the theoretical upper limit based on expensive numerical solvers.

[0067] Computational efficiency: Compared with AD-FWI based on traditional numerical solvers, this method improves the computational speed by about 12 times while maintaining the accuracy of the inverted structure (from about 720 seconds / sample to about 60 seconds / sample). Example 5

[0068] Scalability verification of large-scale 3D forward modeling This embodiment verifies the scalability of the disclosed method in large-scale 3D forward modeling.

[0069] Experimental setup: The experiment was conducted on the Overthrust 3D model with the mesh size expanded to 128×128×128. An axial-AFNO architecture based on axial decomposition was used for modeling.

[0070] Results analysis: High-dimensional adaptability: As shown in Figure 11, compared with the PINO method, the method disclosed in this paper exhibits superior stability in complex 3D media. When the model complexity increases, the prediction results of PINO show severe phase distortion (relative error > 0.86), while the method disclosed in this paper maintains a low relative error of 0.15 ~ 0.19.

[0071] Wavefield characteristics: 3D wavefield slices show that this method accurately recovers the propagation characteristics and fine scattering details of spherical waves in three-dimensional space. This proves that the "axial decomposition + shared basis" architecture adopted in this invention can effectively overcome the memory and computational bottlenecks of 3D modeling.

[0072] While some embodiments of this disclosure have been shown and described, those skilled in the art will understand that modifications may be made to these embodiments without departing from the principles and spirit of this disclosure, which are defined by the claims and their equivalents.

Claims

1. A self-supervised training method for neural operator networks used in solving wave equations, characterized in that, Includes the following steps: Step 1: Define and initialize the neural operator network for processing the probed wave velocity field of matter or object; construct the wave velocity field dataset and determine the source location and frequency to be calculated for each wave velocity field; determine the dynamic iteration step k of BS in physics loss, and the computational parameters required for BS. Step 2: Combine a wave velocity field, a seismic source location, and a frequency obtained through dynamic random sampling into an input sample, and take a set of input samples from the dataset to form a training batch. Step 3: Preprocess the training batch data to make it into a suitable format and input it into the neural operator network; Step 4: Calculate the forward propagation of the neural network and output the predicted wave field; Step 5: Using the predicted wave field as the initial field for BS iteration, perform multiple BS iterations based on the dynamic iteration step number k to output multiple residual fields. Step 6: Take the L2 norm of the k residual fields and sum them to output the physics loss value; Step 7: Determine whether the physics loss has reached the convergence condition. If it has, exit the loop and output the trained neural operator network. If it has not reached the convergence condition, proceed to step 8. Step 8: Call the automatic differentiation function provided by the programming framework to calculate backpropagation and obtain the gradient of the neural operator network parameters; Step 9: Call the optimizer to update the parameters of the neural operator network, and return to step 2 to repeat the process.

2. The self-supervised training method for neural operator networks for solving wave equations as described in claim 1, characterized in that, The specific method for obtaining the frequency through dynamic random sampling is as follows: the training is divided into several stages, and a different frequency sampling range is set for each stage. When the training in a certain stage converges, the next stage is entered.

3. The self-supervised training method for neural operator networks for solving wave equations as described in claim 1, characterized in that, The specific method for multiple BS iterations based on dynamic iteration step k is as follows: the training is divided into several stages, and the BS iteration step k of the physical loss function is dynamically adjusted with each training round. Different k values ​​are set for each stage. When the training in a certain stage converges, the next stage of training begins.

4. The self-supervised training method for neural operator networks for solving wave equations as described in claim 1, characterized in that, The neural operator network is structured as follows: it adopts an axially adaptive Fourier neural operator architecture, specifically including: The input data first passes through an input smoothing projection layer, then through a feature block and position encoding module, before entering the input smoothing projection layer. The core backbone network is composed of layer-axial-AFNO modules; The axial-AFNO module introduces a shared projection substrate, which is shared between different spatial axes in the input data. At the same time, a specific gating mechanism is introduced for each axis. The features are modulated in the frequency domain and then restored by inverse Fourier transform. The smooth projection layer contains stacked convolutional kernels with sizes of [sizes not specified]. 3D convolution.

5. The self-supervised training method for neural operator networks for solving wave equations as described in claim 4, characterized in that, The specific implementation method of introducing a shared projection substrate into the axial-AFNO module is as follows: For the multidimensional features of the input, the network constructs a globally shared, learnable projective weight matrix. ; Specifically, after the input features undergo a one-dimensional real-valued Fast Fourier Transform along a specific spatial axis, they are divided along the channel dimension. Each group contains 1 independent group. There are 1 channel; for the extracted low-frequency component, let its feature within a given group be . , where L is the number of cutoff frequency modes, and the shared basis is linearly reduced and increased in dimensionality through low-rank projection; Separate features into real parts and the virtual part Subsequently, shared matrices were used. , Given a set low-rank size, and k as the current layer number, perform spatial orientation-independent projection calculations: After subsequent gated modulation, a shared matrix is ​​then used. Map the rank dimension back to the feature dimension.

6. The self-supervised training method for neural operator networks for solving wave equations as described in claim 5, characterized in that, The specific implementation of the aforementioned gating mechanism is as follows: A specific gating mechanism is introduced for each axis, specifically by assigning specific, learnable gating parameters and biases to each independent spatial axis in the frequency domain. Furthermore, this gating mechanism introduces a frequency band partitioning strategy. (The truncated portion...) The low-frequency mode is divided into several frequency bands according to frequency level, for the first... The b-th frequency band of each spatial axis is defined with a specific gate vector as follows: After passing through the shared base After projection, the specific gate vector is used to perform element-wise multiplication modulation with the feature: The modulated features are processed by a nonlinear activation function and the above projection and gating process is repeated before entering the next layer; Finally, the frequency domain features processed independently along each axis are restored to the spatial domain through the inverse discrete Fourier transform. When multi-axis feature fusion is used, learnable fusion weight parameters are initialized for each spatial axis. After normalization by the Softmax function, the output wavefields of each axis are weighted and summed to obtain the final wavefield prediction result that balances computational efficiency and anisotropic features.

7. An end-to-end full waveform modeling and updating method, comprising four processing stages: Phase 1: Physics-driven input construction, starting with the velocity model of the input. Random sampling and enhancement are performed, followed by combining the angular frequencies of the random samples. and the location of the epicenter Perform dimensionless processing to generate a normalized velocity model. and normalized source mask ; Phase Two: Parameter-efficient neural prediction. The normalized input first passes through a pre-computation module, which analytically calculates a dimensionless background wave field. To address source singularities, the input data is then fed into a frequency band partitioning expert system. This system comprises multiple parallel expert sub-models, each outputting a predicted value for the scattered wavefield. Finally, the background wavefield is superimposed to obtain the full wavefield prediction. ; Phase 3: CBS-guided physical supervision, constructing an unlabeled physical loss function, and converting the predicted wave field output by the network into a physical loss function. As initial values, a convergent Born series is iterated with a fixed number of steps, and the CBS operator is used. Using the fixed-point property of the wave field, calculate the consistency error between the predicted wave field and the iterated wave field; Phase Four: Velocity inversion based on reflection sensing. The neural operator network parameters trained using the self-supervised training method for wave equation solving as described in any one of claims 1-6 are frozen and used as a differentiable forward modeling engine, introducing the up-wave separation operator. When calculating the backpropagation gradient, the energy of the uplink reflected wave is explicitly separated and weighted, thereby enabling the recovery of the fine underground structure from strong direct wave interference when updating the velocity model.

Images