Seismic parameter inversion method based on RBF (Radial Basis Function) proxy model and multi-algorithm collaborative optimization

Abstract

The invention provides a seismic parameter inversion method based on an RBF proxy model and multi-algorithm collaborative optimization, and relates to the technical field of seismic services, and the method comprises the steps: building an RBF proxy model, carrying out the local real evaluation and proxy verification, and dynamically dividing a global exploration stage and a local refinement stage according to the optimal fitness change rate, smooth transition is carried out by adopting a PSO particle swarm inertia weight linear decline strategy; in the global and local stages, the annealing algorithm inferior solution acceptance rate and the cross entropy sampling distribution covariance are continuously monitored, and dynamic adjustment is carried out. According to the method, rapid prediction of complex stratum parameters and observation waveforms can be achieved, misleading search caused by agent failure is avoided, wide global exploration can be conducted in the early stage through multi-algorithm collaboration and dynamic division of the optimization stage, local refinement can be achieved in the later stage, search activity and diversity are effectively maintained, and the search efficiency is improved. The search efficiency and the inversion reliability are remarkably improved, and meanwhile, the dependence on expensive forward calculation is reduced.

Description

Technical Field

[0001] This invention relates to the field of earthquake service technology, and in particular to a method for earthquake parameter inversion based on the RBF proxy model and multi-algorithm collaborative optimization. Background Technology

[0002] Seismic parameter inversion uses directly observable seismic records to infer the physical parameters of the subsurface medium. This involves observing the propagation characteristics of seismic waves in the subsurface medium to infer the physical parameters of the formation, including P-wave velocity, S-wave velocity, and density. These parameters determine the subsurface structural characteristics, reservoir properties, fault structure, and subsurface medium variations. Since the subsurface cannot be directly observed, core sampling is limited and extremely costly. Therefore, it is necessary to rely on seismic inversion to infer the subsurface structure across the entire region, enabling applications such as seismic exploration, oil and gas reservoir prediction, engineering geological assessment, and seismic hazard analysis.

[0003] For example, Chinese invention patent CN119291774A discloses a method for frequency-varying AVO inversion and parameter prediction of orthogonal fractured reservoirs, including: deriving and reconstructing a fluid-saturated horizontal-vertical orthogonal fractured medium rock physics model based on multiple sets of fracture-pore media Chapman models, realizing fluid-saturated rock physics modeling of media containing two sets of orthogonal fractures; realizing frequency-varying AVO forward modeling of fractured reservoirs based on anisotropic reflectivity algorithm and fluid-saturated horizontal-vertical orthogonal fractured medium rock physics model, and calculating the frequency-varying AVO reflection coefficient and synthetic seismic record of fractured reservoirs under the interface reflection model; driving the frequency-varying AVO forward modeling process, constructing a frequency-varying AVO inversion theory based on a Bayesian framework, utilizing pre-stack seismic AVO data, and introducing simulated annealing-particle swarm optimization algorithm in the inversion process to realize frequency-varying AVO inversion based on a Bayesian inversion framework and predict reservoir parameters.

[0004] For example, the Chinese invention patent with announcement number CN106842310B discloses a pre-stack seismic four-parameter synchronous inversion method, which includes: based on the viscoelastic medium theory, the method establishes a direct relationship between the seismic wave reflection coefficient and the formation P-wave velocity, S-wave, density and formation absorption parameter Q, and realizes the synchronous inversion of the four parameters of formation P-wave velocity, S-wave, density and formation absorption parameter Q through a stable pre-stack seismic inversion algorithm.

[0005] Traditional seismic parameter inversion methods typically rely on numerical simulation techniques such as least squares, regularized inversion, or full waveform inversion. With the development of surrogate models, machine learning, and heuristic optimization algorithms, seismic inversion methods have begun to introduce methods that construct surrogate models based on observed seismic waveform data. This maps formation parameters to corresponding waveform features and combines particle swarm optimization and annealing algorithms to search for optimal solutions in the parameter space. During the algorithm iteration process, the algorithm feeds back the optimization results and updates the surrogate model, enabling the model and the search algorithm to form a closed-loop collaboration.

[0006] However, when combining particle swarm optimization (PSO) with simulated annealing, it's possible that the samples generated during the optimization process of both algorithms in the high-dimensional parameter space are concentrated near a few local optima. This lack of sample diversity makes it easy for the algorithms to overlook potential global optimum regions. To increase sample diversity and broaden the search area, cross-entropy is introduced to update the probability distribution of diverse samples. However, because the annealing temperature of the annealing algorithm often lags behind the information sharing between PSO and cross-entropy due to differences in update mechanisms and step sizes, the rate of annealing temperature decrease becomes uncontrolled. The probability of the annealing algorithm accepting inferior solutions rapidly decreases, retaining only those solutions that currently appear better. It loses the ability to escape the current algorithm's detection area by accidentally accepting inferior solutions. Simultaneously, when updating the probability distribution using these limited better solutions, the cross-entropy algorithm may mistakenly believe that optimal solutions have already appeared due to the lack of sample diversity, further narrowing the covariance range of the distribution and forming an overly concentrated sampling distribution. The search space is compressed, and the annealing algorithm, due to its excessively low annealing temperature, cannot break this local concentration. The entire optimization process is locked in a local optimum region. Summary of the Invention

[0007] Therefore, this invention provides a seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization, which can provide reliable optimization results for seismic exploration and stratigraphic parameter inversion.

[0008] This invention provides a seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization. The method includes:

[0009] Raw seismic observation data is collected, standardized, and then an RBF surrogate model is constructed for inversion and iterative optimization.

[0010] During the inversion iteration process, the timestamps of the shared states between algorithms are collected in real time, proxy verification is performed, the adjustment requirements of the RBF proxy model are determined, and the optimization stage is dynamically divided.

[0011] During the optimization phase, which is a global exploration phase, the annealing temperature adjustment strategy and the cross-entropy distribution adjustment strategy are determined.

[0012] When the optimization phase is a local refinement phase, a multi-algorithm collaborative adjustment strategy is determined and implemented.

[0013] After receiving the inversion execution signal, a closed-loop optimization inversion is performed to obtain the globally optimal inversion result.

[0014] The beneficial effects of the technical solutions provided in the embodiments of the present invention include at least the following:

[0015] 1. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization provided by this invention achieves rapid prediction of complex stratigraphic parameters and observed waveforms by constructing an RBF surrogate model, avoiding misleading searches caused by surrogate failure. Through multi-algorithm collaboration and dynamic division of optimization stages, extensive global exploration can be carried out in the early stage to quickly locate potential superior areas, and local refinement can be achieved in the later stage to improve convergence speed and accuracy. Adaptive adjustment of annealing temperature, cross-entropy sampling distribution expansion, short-term tabu mechanism, and random sample injection effectively maintain search activity and diversity, avoiding the algorithm from getting stuck in local lock-in or premature convergence. Through rapid prediction, accuracy control, multi-algorithm collaboration, and adaptive adjustment mechanism, efficient, stable, and accurate joint inversion of stratigraphic parameters is achieved, significantly improving search efficiency and inversion reliability, while reducing dependence on expensive forward modeling calculations.

[0016] 2. This invention dynamically controls the probability of the algorithm accepting inferior solutions by adjusting the annealing temperature, thereby maintaining sufficient exploration capability in the global search phase and avoiding premature convergence to local optima. In the local refinement phase, the search intensity is finely adjusted by adaptively lowering or raising the temperature, enabling particles to converge fully in the optimal region. At the same time, when the search stagnates or becomes locally locked, the activity is restored by rolling back to the stable temperature and superimposing short-term perturbations. This ensures that the algorithm can both extensively explore the parameter space and efficiently converge to the optimal solution at different stages, improving the stability and reliability of the optimization process.

[0017] 3. This invention adaptively controls the distribution range and concentration of samples in the search space by adjusting the cross-entropy sampling distribution. It expands the sampling coverage and increases search diversity by expanding the covariance matrix, preventing local convergence caused by excessive concentration of candidate solutions. At the same time, by adjusting the elite ratio to balance the dependence on excellent samples, the algorithm can fully explore the global potential excellent region and focus on high fitness regions for detailed search during the local refinement stage. This improves optimization efficiency, enhances global search capability, and effectively reduces the risk of getting trapped in local optima.

[0018] 4. This invention promptly identifies regions that are stuck or continuously failing during the search process through local lockout adjustment, and restores search activity and diversity by increasing the annealing temperature, freezing the descent rounds, injecting new random samples, and enabling short-term tabu mechanisms. This prevents the algorithm from staying in local optima for a long time, promotes the search to migrate to unexplored regions, and achieves effective optimization on a global scale. This not only improves the robustness and adaptability of the algorithm, but also ensures that potential excellent solutions can be continuously discovered during the local refinement stage, thereby improving the overall inversion accuracy and reliability. Attached Figure Description

[0019] Figure 1This is a flowchart of the seismic parameter inversion method based on RBF proxy model and multi-algorithm collaborative optimization provided in this embodiment of the invention;

[0020] Figure 2 This is a step-by-step strategy diagram of the seismic parameter inversion method based on the RBF proxy model and multi-algorithm collaborative optimization provided in this embodiment of the invention;

[0021] Figure 3 This is a diagram illustrating the RBF proxy model construction strategy for the seismic parameter inversion method based on the RBF proxy model and multi-algorithm collaborative optimization provided in this embodiment of the invention.

[0022] Figure 4 This is a diagram illustrating the multi-algorithm collaborative adjustment strategy of the seismic parameter inversion method based on the RBF proxy model and multi-algorithm collaborative optimization provided in this embodiment of the invention. Detailed Implementation

[0023] Particle Swarm Optimization (PSO) achieves rapid global exploration and solution space convergence through collaborative search, efficiently locating potential optimal solution regions in high-dimensional complex models. Annealing (SA) simulates the physical annealing process, accepting inferior solutions with controlled probability to escape local optima traps and enhance the robustness and stability of global search. Cross-entropy sampling (CE) continuously updates probability distribution parameters, guiding samples towards high-quality solution sets, ensuring the search distribution gradually focuses on the optimal region. The combination of these three algorithms forms an optimization system integrating global exploration, random perturbation, and distributed convergence, achieving efficient search, robust jumps, and accurate convergence in the inversion process, improving the accuracy and reliability of model parameter estimation.

[0024] To make the objectives, technical solutions, and advantages of this application clearer, the embodiments of this application will be described in further detail below with reference to the accompanying drawings.

[0025] In this application, the terms "first," "second," and "third," etc., are used to distinguish identical or similar items with essentially the same function. It should be understood that there is no logical or temporal dependency between "first," "second," and "nth," nor does it limit the quantity or execution order. It should also be understood that although the following description uses the terms "first," "second," etc., to describe various elements, these elements should not be limited by the terms. These terms are merely used to distinguish one element from another; for example, the first local refinement criterion and the second local refinement criterion are only used to distinguish sample data.

[0026] It should also be understood that, in the various embodiments of this application, the sequence number of each process does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0027] It should also be understood that determining B based on A does not mean determining B solely based on A; it is also possible to determine B based on A and / or other information.

[0028] It should also be understood that the term “comprising” (also referred to as “includes”, “including”, “comprises” and / or “comprising”) as used in this specification specifies the presence of the stated features, integers, steps, operations, elements, and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof.

[0029] It should also be understood that the term "if" can be interpreted as meaning "when" or "upon" or "in response to determination" or "in response to detection." Similarly, depending on the context, the phrases "if determination..." or "if detection [the stated condition or event]" can be interpreted as meaning "when determination..." or "in response to determination..." or "when detection [the stated condition or event]" or "in response to detection [the stated condition or event]."

[0030] like Figure 1 The flowchart shown is for a seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization. The processing flow of this method can include the following steps: constructing a training set using stratigraphic parameter vectors and observed waveform data, establishing an RBF surrogate model, and performing local real-world evaluation and surrogate verification by real-time monitoring of shared state timestamps and optimal fitness change rates through each module. Simultaneously, the global exploration and local refinement stages are dynamically divided according to the optimal fitness change rate, and a linear decreasing strategy of PSO particle swarm inertia weight is adopted to achieve a smooth transition from exploration to refinement. In the global and local stages, the acceptance rate of inferior solutions and the cross-entropy sampling distribution covariance of the annealing algorithm are continuously monitored and dynamically adjusted. When local deadlock occurs, an emergency recovery mechanism is triggered to inject random samples, enable short-term taboos and rapid temperature control responses to maintain search diversity and activity.

[0031] like Figure 2The diagram illustrates the steps of a seismic parameter inversion method based on an RBF surrogate model and multi-algorithm collaborative optimization. These steps include: establishing a surrogate model using a cubic spline RBF kernel function; obtaining weight parameters by solving a system of linear equations; real-time monitoring of shared state timestamps and the rate of change of optimal fitness in each module; performing local real-world assessment and surrogate verification; generating seismic waveforms and calculating the mean square error through forward modeling; adjusting the RBF kernel width or refreshing the local state accordingly; dynamically dividing the process into global exploration and local refinement stages based on the rate of change of optimal fitness; and utilizing a linearly decreasing PSO particle swarm inertial weight strategy to achieve a smooth transition from wide-area search to local convergence. During both the global and local stages, continuously monitoring the acceptance rate of inferior solutions in the annealing algorithm and the covariance of the cross-entropy sampling distribution, and implementing annealing temperature adjustment, cross-entropy distribution adjustment, and local lockout adjustment strategies. After strategy execution, the surrogate model is retrained with a delay to adapt to the new sample distribution, ensuring that the prediction accuracy is consistent with the search direction, thereby achieving joint inversion of formation parameters and acquisition of multi-objective optimal solutions.

[0032] like Figure 3 The diagram illustrates the RBF surrogate model construction strategy for seismic parameter inversion based on the RBF surrogate model and multi-algorithm collaborative optimization. This strategy includes: standardizing the original seismic observation data and determining a reasonable range of inversion parameters based on prior geological knowledge; then using Latin hypercube sampling technology to generate representative and uniformly distributed initial sample points in the parameter space; performing forward modeling on each sample point and evaluating its fit with the observation data; selecting the optimally fitted sample point as the center point of the radial basis function RBF surrogate model; constructing the RBF model using cubic spline basis functions; and determining its optimal weight parameters by solving a system of linear equations.

[0033] Raw seismic observation data is collected, standardized, and then an RBF surrogate model is constructed for inversion and iterative optimization.

[0034] Furthermore, an RBF proxy model is constructed, and the specific analysis method is as follows:

[0035] The raw seismic observation data is standardized, and the reasonable range of inversion parameters is determined based on prior geological knowledge. The Latin hypercube sampling technique is used to generate initial sample points that are representative and uniformly cover the entire parameter space.

[0036] It is important to explain that raw seismic observation data mainly includes time-series ground vibration signals acquired by accelerometers, velocity meters, or displacement sensors. These signals record metadata such as ground acceleration, velocity, or displacement changing over time. Different types of seismic waves can be identified from this data: the first arriving P-wave (P-wave) and the subsequent arriving S-wave (S-wave), reflecting the propagation characteristics of seismic energy. Geological prior knowledge refers to the known or empirical information about the subsurface structure, lithological distribution, sedimentary environment, and range of physical parameters of the study area, based on existing geological, geophysical, and geochemical data, prior to seismic inversion or stratigraphic modeling. This includes stratigraphic division, fault distribution, lithological types, well logging data, velocity models, and density ranges. This knowledge provides parameter constraints and a reasonable initial model for inversion, helping to narrow the search space, improve the geological rationality and stability of the model, and avoid inversion results that are physically impossible or inconsistent with geological laws. Therefore, it plays a guiding and corrective role in the algorithm optimization process, making the inversion results more consistent with actual geological characteristics.

[0037] It should be noted that the raw seismic observation data undergoes standardization processing. Based on prior geological knowledge, reasonable ranges for the inversion parameters are determined, including upper and lower limits for P-wave velocity, S-wave velocity, and density. Latin hypercube sampling is used to generate representative initial sample points within the parameter space. These sample points uniformly cover the entire parameter space. The inversion parameter space is defined as follows: the search boundaries for P-wave velocity Vp, S-wave velocity Vs, and density ρ are defined. An initial sample set S = {S1, S2, ... S} is generated using Latin hypercube sampling. M}, where each sample S i Representing a set of formation parameters (Vp, Vs, and ρ), the sample size M is set, and the true fitness value of the initial sample is evaluated by means of waveform fitting error.

[0038] Forward modeling is performed on each initial sample point to evaluate its fit with the observed data.

[0039] The best-fitting sample point is selected from the initial sample points as the center point of the radial basis function model. The RBF surrogate model is constructed using cubic spline basis functions. The optimal weight parameters of the RBF surrogate model are determined by solving the linear equation system.

[0040] It should be explained that the fitting effect is achieved by comparing the difference between the simulated seismic waveform and the actual observed seismic waveform. In practice, after inputting the stratigraphic parameter vector into the forward model to generate the simulated seismic waveform, the mean square error between the simulated waveform and the observed waveform is calculated. The smaller the mean square error, the higher the degree of fit between the two, and the closer the simulation result is to the real observation data. The fitting index of all initial sample points is sorted, and the sample with the smallest mean square error is selected as the sample with the best fitting effect, which represents the solution that best matches the current stratigraphic parameters and the observation data. This solution is used for the selection of the center point of the surrogate model or parameter update in the future.

[0041] It needs to be explained that the forward model calculates the seismic response based on known geological physical parameters (such as P-wave velocity, S-wave velocity, and density), and its output is usually a ground vibration signal in the form of a time series. The forward model and the inversion model are essentially corresponding to each other. The forward model predicts the waveform from the geological parameters, while the inversion model estimates the geological parameters from the observed waveform.

[0042] The radial basis function model is defined as follows:

[0043] ,

[0044] f(x) represents the local response function of the RBF surrogate model, Wi is the inertia weight, and ϕ(r) = r 3 Let r be the cubic spline basis function, ci be the Euclidean distance, ci be the RBF center point, and a sample be selected from the initial sample as the center. P(x) is a linear polynomial term used to capture the global trend, satisfying the following condition: Let i represent the i-th sample, i = 1, 2, 3, ..., M.

[0045] It should be explained that the RBF surrogate model constructed using cubic spline basis functions first selects the input and output datasets of initial sample points and normalizes the input variables to eliminate the influence of dimensions. Then, several nodes in the input space are selected as the support centers of the spline basis, and the cubic spline basis functions are used as the kernel functions to characterize the nonlinear mapping relationship between the input variables and the output response. The model establishes the coefficient matrix by minimizing the sample prediction error, and solves for the weight coefficients of each spline basis using regularization or least squares methods, achieving a balance between global smoothness and local fine-fitting. The resulting RBF surrogate model exhibits good continuity and differentiability, and can achieve high-precision approximation and stable generalization in complex nonlinear problems.

[0046] When constructing an RBF surrogate model, the optimal weight parameters for each radial basis function need to be determined by solving a system of linear equations. Specifically, the following steps are taken: First, a radial basis function matrix is ​​constructed using the training samples as input, where each element is obtained by transforming the Euclidean distance between samples using a kernel function. Then, this matrix is ​​mapped to the output vector to be fitted, forming a system of linear equations Φw=y, where Φ is the radial basis function matrix, w is the weight coefficient vector to be determined, and y is the true output value. The optimal solution can be obtained using the least squares method, thus obtaining the best weight parameters that balance fitting accuracy and model stability, achieving high-precision approximation of complex nonlinear relationships using the RBF surrogate model.

[0047] During the inversion iteration process, the timestamps of the shared states between algorithms are collected in real time, proxy verification is performed, the adjustment requirements of the RBF proxy model are determined, and the optimization stage is dynamically divided.

[0048] Furthermore, perform proxy verification to determine the adjustment requirements of the RBF proxy model. The specific analysis method is as follows:

[0049] An initial sample set is constructed for each initial sample point, and then iterative optimization is performed using multiple algorithms. These algorithms include: PSO algorithm, simulated annealing algorithm, and cross-entropy algorithm.

[0050] It should be noted that in the inversion iterative optimization process, each initial sample point is first aggregated to form an initial sample set. These samples uniformly cover the entire parameter space and represent possible combinations of formation parameters. The initial sample set is not only used to train the initial RBF surrogate model but also serves as the starting basis for multi-algorithm iterative optimization. In the subsequent optimization process, the PSO algorithm, simulated annealing algorithm, and cross-entropy algorithm use the initial sample set as a starting point to iteratively generate new candidate formation parameter vectors. They also combine the RBF surrogate model to quickly predict the objective function value to evaluate the quality of the samples, thereby dynamically updating the sampling distribution during the global search and local refinement stages, gradually approaching the optimal inversion solution.

[0051] The timestamps of shared states between algorithms are collected in real time. When it is detected that the timestamp read by any algorithm is inconsistent with the current timestamp, the cross-algorithm data interaction of the corresponding algorithm is temporarily suspended, and proxy verification is performed first. The specific analysis process is as follows:

[0052] The current iteration's stratigraphic parameter vector is obtained as a sample and input into the forward modeling module for simulation to obtain the simulated seismic waveform.

[0053] It should be noted that during the inversion iterative optimization process, each algorithm update generates a set of candidate stratigraphic parameter vectors. These vectors can be understood as samples generated in the current iteration. Each sample is input into the forward modeling module for simulation. The corresponding seismic waveform response is calculated using the seismic wave propagation model. The generated simulated waveform reflects the response of the set of stratigraphic parameters to the observed wave under the theoretical model. Subsequently, it can be compared with the actual observed seismic waveform, and the quality of the sample can be evaluated by calculating the error or goodness of fit.

[0054] It is important to explain that the seismic wave propagation model is used to describe the propagation characteristics of seismic waves in subsurface media and is the theoretical foundation for geophysical inversion and seismic imaging. This model treats the strata as a heterogeneous medium composed of different elastic parameters (such as density, P-wave velocity, and S-wave velocity). By solving the wave equation, it simulates the propagation, reflection, refraction, and scattering processes of seismic waves in the medium. The model can adopt either the acoustic wave equation or the elastic wave equation to describe the propagation characteristics of P-waves and S-waves, respectively. Through numerical calculations or analytical approximations, the seismic wave propagation model can reflect the influence of subsurface structures on the wave field, providing a physical basis for stratum parameter inversion and geological structural interpretation.

[0055] The simulated seismic waveform is compared with the observed seismic waveform, and the mean square error is calculated.

[0056] It should be noted that the simulated seismic waveform and the observed seismic waveform are aligned in time or sampling points to ensure that the corresponding sampling points match one by one. Then, for each time point, the difference between the simulated value and the observed value is calculated, and the difference is squared to eliminate the sign effect. The sum of the squared differences of all time points and divided by the total number of sampling points gives the mean square error of the sample.

[0057] Extract the preset mean square error threshold from the database.

[0058] It's important to explain that the preset mean squared error (MSE) threshold is typically determined based on a combination of historical sample statistics and model validation results. First, a large number of samples with known inputs and actual outputs are extracted from the database. An RBF surrogate model or other predictive model is initially trained, and the MSE distribution of each sample is calculated. Then, statistical analysis is performed on the error distribution, such as calculating its mean, standard deviation, or quantiles, to reflect the typical error range of the model under different complexity or noise levels. Combining empirical knowledge or inversion accuracy requirements, a MSE threshold is set based on this, for example, taking the mean plus a certain number of standard deviations, or taking a specific upper quantile of the distribution, as a control standard for model accuracy. This threshold is stored in the database and used as a basis for subsequent model evaluation and adaptive adjustment.

[0059] If the mean squared error is greater than or equal to the mean squared error threshold, the adjustment requirement of the RBF proxy model is recorded as demand adjustment. The RBF kernel width is adjusted to correct the proxy model by the deviation between the mean squared error and the mean squared error threshold, and the local state cache is refreshed.

[0060] It should be noted that if the mean squared error is greater than or equal to the mean squared error threshold, it means that the surrogate model's prediction accuracy for this parameter combination is insufficient and there is a large error, requiring model correction. In this case, the RBF surrogate model adjustment requirement is recorded as a requirement adjustment, indicating that the current surrogate model needs to be updated. The mean squared error deviation value is obtained by subtracting the mean squared error threshold from the mean squared error to determine the RBF kernel width adjustment value to correct the surrogate model. After the kernel width adjustment is completed, the updated model parameters are saved and the local state cache is refreshed to ensure that subsequent iterations use the latest surrogate model and state information, thereby ensuring that the optimization search direction is consistent with the surrogate prediction accuracy.

[0061] The RBF core width adjustment value is added to the current RBF core width in the system log to obtain the updated RBF core width value, and then the updated RBF core width value is applied to the next iteration.

[0062] In this embodiment, a mapping rule between the mean squared error deviation value and the RBF kernel width adjustment value is established in advance and stored in a configuration file in the form of a mapping table or function model for unified invocation and automated execution during multi-algorithm optimization. This mapping rule reflects the adaptive adjustment logic of the kernel width of the surrogate model under different error levels: when the deviation value is large, it indicates that the prediction accuracy of the surrogate model in the current parameter region is insufficient, and the kernel width needs to be increased to expand the response range of the radial basis function, thereby improving the local fitting ability; when the deviation value is small or close to the threshold, the kernel width can remain unchanged or be slightly reduced. In actual operation, the optimization algorithm first calculates the mean squared error between the simulated waveform and the observed waveform of the current sample and compares it with the preset mean squared error threshold to obtain the mean squared error deviation value. This mean squared error deviation value is then input into the mapping table or function model to find or calculate the corresponding kernel width adjustment amount, which is then applied to the kernel function parameters of the RBF surrogate model to correct the model's prediction ability in that parameter region.

[0063] It should be explained that the RBF kernel width is a parameter that controls the range and influence of the radial basis function response. Adjusting the RBF kernel width can make the radial basis function prediction more accurate in the current parameter region.

[0064] If the mean squared error is less than the mean squared error threshold, the adjustment requirement of the RBF proxy model is marked as no adjustment is needed, the use of this sample for cross-algorithm updates is temporarily suspended, and a local state refresh is triggered.

[0065] It should be noted that if the mean squared error is less than the mean squared error threshold, it means that the prediction accuracy of the surrogate model in the current parameter region has met the requirements, and there is no need to further adjust the kernel width or correct the model. Therefore, the adjustment requirement of the RBF surrogate model is recorded as no adjustment is needed. This sample is temporarily suspended for cross-algorithm data updates to avoid unnecessary interference to the global search distribution. At the same time, a local state refresh is triggered to update the local cache information and ensure that the surrogate model parameters and state used in subsequent iterations are up-to-date.

[0066] It should be noted that local state refresh refers to updating the algorithm's internal state variables and control parameters in real time based on the latest calculation results during optimization or inversion. This maintains dynamic consistency and information freshness in the search process. First, key metrics generated in the current iteration, such as the objective function value, search particle distribution, acceptance rate, temperature, or sampling variance, are collected and compared with the previous state record. If significant changes or performance degradation are detected, a state refresh operation is triggered. The refresh process includes: clearing or resetting part of the historical cache, updating the current optimal solution and its corresponding parameters, synchronizing the latest prediction results of the proxy model, and reinitializing some search individuals or adjusting the weight distribution if necessary.

[0067] Furthermore, a dynamic partitioning and optimization phase is conducted, with the specific analysis method as follows:

[0068] Calculate the rate of change of the optimal fitness in each iteration.

[0069] It should be noted that the fitness value of a sample is obtained by taking the reciprocal of its mean squared error. The fitness values ​​of all candidate samples in the current iteration are recorded, and the optimal fitness (the sample with the highest fitness) is determined. This is compared with the optimal fitness from the previous iteration, and the change in optimal fitness is calculated as the current optimal fitness minus the previous optimal fitness. This change is then divided by the iteration step size or the corresponding time interval to obtain the rate of change of optimal fitness. This rate reflects the speed at which the optimal value of the objective function improves per unit iteration or unit time.

[0070] Extract the preset optimal fitness change rate threshold from the database.

[0071] When the rate of change of the optimal fitness is greater than or equal to the rate of change threshold, the optimization phase is divided into the early global exploration phase.

[0072] It should be noted that when the rate of change of the optimal fitness is greater than or equal to the rate of change threshold, it means that the optimal fitness has been significantly improved in the current iteration, the search is still in a state of rapid improvement, indicating that the algorithm has not yet converged and is extensively exploring the parameter space. Based on this judgment, the current optimization stage can be divided into the early global exploration stage. The focus of this stage is to expand the search range, quickly locate potential good areas, and maintain search diversity in order to avoid getting trapped in local optima too early.

[0073] When the rate of change of the optimal fitness is less than the rate of change threshold, the optimization phase is divided into the late local refinement phase.

[0074] It should be noted that when the rate of change of the optimal fitness is less than the rate of change threshold, it means that the improvement of the optimal fitness in the current iteration is limited and the search speed is significantly slowed down. This indicates that the algorithm is close to the excellent region or the local convergence state. Based on this judgment, the current optimization stage can be divided into the late local refinement stage. The focus of this stage is to conduct a fine search and continuous convergence of the discovered potential excellent regions to improve the accuracy of candidate solutions.

[0075] Furthermore, adaptive adjustment of the PSO algorithm parameters is performed, and the specific analysis method is as follows:

[0076] Extract the preset initial inertia weight, final inertia weight, and maximum number of iterations from the database.

[0077] It should be noted that during the iterative optimization process, in order to adapt to the local refinement requirements, the inertia weight of the PSO particle swarm optimization algorithm is usually gradually reduced according to a linear decreasing strategy, so that the particles gather in the local good region, improve the convergence ability, and at the same time ensure the stability and efficiency of the search process. The inertia weight controls the degree of inertia of the particle at the current position, that is, the ability of the particle to maintain its original velocity in the search space.

[0078] It should be noted that the database pre-sets the initial inertia weight, final inertia weight, and maximum number of iterations for the particle swarm optimization algorithm. This pre-setting logic is designed based on the algorithm's search requirements at different stages: a larger initial inertia weight enhances the global exploration capability of the particle swarm, giving particles higher inertia in the early stages of the search to cover a wider solution space; a smaller final inertia weight improves local convergence accuracy in the later stages of iteration, allowing particles to gradually focus on high-quality solution regions; and the maximum number of iterations is determined based on problem complexity and computational resource constraints to balance computational overhead and convergence performance. This logic ensures the algorithm maintains sufficient global exploration activity in the early stages and achieves stable convergence in the later stages, thereby improving overall inversion and optimization performance.

[0079] Based on the ratio of the current iteration count to the maximum iteration count in the system log, the current inertia weight is determined to be the ratio between the initial inertia weight and the final inertia weight, thus determining the current inertia weight of the PSO particle swarm.

[0080] It should be noted that, Where Wmax is the initial inertia weight, Wmin is the final inertia weight, Tmax is the maximum number of iterations, t is the current number of iterations, and W is the current inertia weight of the PSO particle swarm.

[0081] During the optimization phase, which is a global exploration phase, the annealing temperature adjustment strategy and the cross-entropy distribution adjustment strategy are determined.

[0082] Furthermore, the annealing temperature adjustment strategy and the cross-entropy distribution adjustment strategy were determined, and the specific analysis methods are as follows:

[0083] In the early global exploration phase, the acceptance rate of inferior solutions of the annealing algorithm is collected and recorded as the first inferior solution acceptance rate of the annealing algorithm. At the same time, the cross-entropy sampling distribution covariance is collected and recorded as the first cross-entropy sampling distribution covariance.

[0084] It should be explained that the acceptance rate of inferior solutions in the annealing algorithm is obtained by statistically analyzing the proportion of candidate solutions whose fitness is lower than the historical best solution in the current iteration to the total number of accepted solutions, while the covariance of the cross-entropy sampling distribution is obtained by calculating the outer product of the deviations of the parameter vector and the mean vector of the elite sample set and taking the average.

[0085] Extract the pre-defined inferior solution acceptance rate threshold and covariance threshold from the database.

[0086] It should be noted that the database pre-sets a suboptimal solution acceptance rate threshold for the annealing algorithm and a covariance threshold for the cross-entropy sampling algorithm. The pre-setting logic is designed based on the principle of balancing global search activity and distribution stability: the suboptimal solution acceptance rate threshold measures the algorithm's tolerance for suboptimal solutions during the exploration phase; a higher threshold enhances the ability to escape local optima and improves global search diversity. Conversely, a lower threshold helps strengthen solution selection and improve result accuracy during the convergence phase. The covariance threshold determines the convergence degree of the sampling distribution. When the distribution covariance is less than the threshold, it means the samples have concentrated within a small range, and the algorithm can enter the local refinement phase. If the covariance is large, a high search diffusion is maintained to prevent premature convergence. Through this logic, adaptive switching between the search and convergence phases can be achieved, maintaining the stability and flexibility of the optimization process.

[0087] If the acceptance rate of the first inferior solution in the annealing algorithm is lower than the threshold for the acceptance rate of inferior solutions, it is recorded as the first condition for global adjustment.

[0088] It should be noted that if the acceptance rate of the first inferior solution in the annealing algorithm is lower than the inferior solution acceptance rate threshold, it means that the algorithm is currently accepting a poor solution with too low a probability, the search activity is insufficient, the global exploration ability may be limited, and there is a risk of premature convergence or local stagnation. This state is marked as the first condition for global adjustment, which serves as the basis for judging the triggering of the global adjustment strategy, thereby ensuring the diversity and exploration ability of the global search phase.

[0089] If the covariance of the first sampling distribution of cross-entropy is less than or equal to the covariance threshold, it is denoted as the second condition for global adjustment.

[0090] It should be noted that if the covariance of the first sampling distribution of cross-entropy is less than or equal to the covariance threshold, it means that the coverage of the current sampling distribution in the search space is too small, the generated candidate samples are too concentrated, the search diversity is insufficient, and the algorithm is prone to premature convergence or local stagnation. The breadth of the sampling distribution is insufficient, and it cannot fully explore potential good areas. It is necessary to increase the spatial coverage of the samples to improve the global search capability.

[0091] If there is a first condition for global adjustment and there is no second condition for global adjustment, then the adjustment strategy in the global exploration phase is recorded as the annealing temperature adjustment strategy.

[0092] It should be noted that if the first condition for global adjustment exists and the second condition for global adjustment does not exist, it means that the main problem of the current search is that the annealing algorithm has a low probability of accepting inferior solutions, while the sample distribution is still sufficiently broad and has not yet shown a situation of overly concentrated distribution or insufficient coverage. The algorithm only needs to adjust the annealing temperature to improve the acceptance rate of inferior solutions and search activity, without adjusting the cross-entropy sampling distribution, thereby maintaining the breadth and diversity of global exploration. Therefore, the adjustment strategy in the global exploration phase is denoted as the annealing temperature adjustment strategy.

[0093] If there is a second global adjustment condition but no first global adjustment condition, then the global exploration phase adjustment strategy is denoted as the cross-entropy distribution adjustment strategy.

[0094] It should be noted that if there is a second global adjustment condition but no first global adjustment condition, it means that the main problem of the current search is insufficient sampling distribution range and the candidate samples are too concentrated, which may lead to insufficient coverage of the exploration space. The annealing algorithm has a normal ability to accept inferior solutions and there is no need to adjust the annealing temperature. In this case, it is necessary to adjust the cross-entropy sampling distribution to increase the search diversity and coverage, while maintaining the activity of the annealing search. Therefore, the adjustment strategy in the global exploration phase is denoted as the cross-entropy distribution adjustment strategy.

[0095] If neither of the two conditions exists, the global exploration phase adjustment strategy is recorded as continuing iterative optimization, and an inversion execution signal is generated.

[0096] It should be noted that the absence of both conditions indicates that the search process is in an ideal state: the algorithm's exploration activity and sample distribution range meet the requirements, and no additional adjustments are needed. The global exploration stage adjustment strategy can be recorded as continuing iterative optimization, keeping the current search strategy unchanged, and sending an inversion execution signal to the optimization execution module, instructing the algorithm to continue generating candidate formation parameter vectors according to the existing parameters and perform forward simulation, thereby realizing continuous inversion calculation and global search.

[0097] If both conditions are met, then the annealing temperature adjustment strategy and the cross-entropy distribution adjustment strategy are executed simultaneously.

[0098] Furthermore, an annealing temperature adjustment strategy is implemented, and the specific analysis method is as follows:

[0099] The inferior solution acceptance rate deviation is obtained by subtracting the inferior solution acceptance rate threshold from the first inferior solution acceptance rate of the annealing algorithm. The reduction value of the annealing temperature decrease rate is determined based on the inferior solution acceptance rate deviation value. The inferior solution acceptance rate of the annealing algorithm is statistically analyzed in real time within a sliding window. The temperature decay coefficient is corrected based on the reduction value of the annealing temperature decrease rate.

[0100] It should be noted that the temperature decay coefficient is added to the decrease in the annealing temperature rate to obtain the updated temperature decay coefficient value, which is then applied to the next iteration. The corrected temperature decay coefficient continues to affect the annealing temperature update in subsequent iterations, making the temperature change more stable and helping the algorithm achieve a balance between global exploration and local refinement.

[0101] It should be noted that in this embodiment, a mapping table or function model between the inferior solution acceptance rate deviation and the annealing temperature decrease rate is pre-established and stored in a configuration file for unified algorithm invocation. The preset logic of this mapping table is as follows: when the inferior solution acceptance rate deviation is large, a larger decrease rate of temperature decrease corresponds to significantly slowing down the annealing temperature decay, thereby increasing the probability of accepting inferior solutions and enhancing global search activity; when the deviation is small or close to the target value, a smaller or zero decrease rate of temperature decrease corresponds to maintain the original rhythm of temperature decay and ensure the stability of the search process. In actual operation, the optimization algorithm first calculates the inferior solution acceptance rate of the current iterative annealing algorithm minus the inferior solution acceptance rate threshold to obtain the inferior solution acceptance rate deviation; then, this deviation value is input into the mapping table or function model to find or calculate the corresponding decrease rate of temperature decrease, and applied to the annealing temperature update, thereby realizing adaptive control of temperature adjustment, enabling the algorithm to dynamically adjust the global exploration intensity under different search states.

[0102] If the acceptance rate of inferior solutions in the adjusted annealing algorithm is still lower than the threshold for the acceptance rate of inferior solutions, the algorithm is rolled back to the previous stable temperature and the superimposed short-term disturbance factor is determined based on the quadratic deviation between the acceptance rate of inferior solutions in the adjusted annealing algorithm and the threshold for the acceptance rate of inferior solutions, in order to restore search activity.

[0103] It should be noted that if the acceptance rate of inferior solutions in the adjusted annealing algorithm is still lower than the threshold, it indicates that the current search activity is still insufficient, and the algorithm is still unable to effectively accept candidate solutions that are worse than the current solution. There is a risk of getting trapped in local optima or the search stagnates. In order to avoid getting trapped in local optima, the annealing temperature is rolled back to the previous stable value to ensure that the temperature is within a reliable range. The short-term perturbation factor is determined by subtracting the threshold from the adjusted annealing algorithm's inferior solution acceptance rate and adding it directly to the previous annealing temperature. This allows the next iteration to temporarily increase the probability of accepting inferior solutions and activate the global search capability.

[0104] In this embodiment, when the acceptance rate of suboptimal solutions in the adjusted annealing algorithm is still lower than the threshold for the acceptance rate of suboptimal solutions in the annealing algorithm, a secondary suboptimal solution acceptance rate deviation value is calculated, which is the difference between the adjusted suboptimal solution acceptance rate and the threshold for the acceptance rate of suboptimal solutions. Typically, a squared operation is used to amplify the impact of the deviation on the perturbation. Based on historical optimization experiments or simulation analysis, the degree of insufficient search activity corresponding to different deviation values ​​and the required perturbation amplitude are statistically analyzed. The secondary deviation value is divided into several intervals, and a corresponding perturbation factor value is defined for each interval. The larger the deviation, the larger the corresponding perturbation factor, so as to ensure that the search activity can be effectively restored. This perturbation factor is used to add a moderate temperature increase or random fluctuation on the basis of rolling back to the previous stable temperature, so as to temporarily increase the probability of accepting suboptimal solutions, thereby activating search activity, expanding the exploration range, and avoiding the algorithm from getting stuck in a local stagnation.

[0105] The specific analysis method for cross-entropy distribution adjustment strategy is as follows:

[0106] Collect the covariance of the first sampling distribution of cross-entropy and extract the covariance threshold from the database.

[0107] If the covariance of the first sampling distribution of cross-entropy is less than or equal to the covariance threshold, then the cross-entropy distribution adjustment strategy is executed, which expands the covariance matrix and reduces the elite ratio. The reduction range of the elite ratio is determined by combining the deviation between the first sampling distribution covariance of cross-entropy and the covariance threshold with the preset maximum allowable reduction range of the elite ratio, thereby obtaining the updated elite ratio value, which is used in the next iteration.

[0108] It's important to explain that the elite ratio, in optimization algorithms based on sample selection or distribution updates, is a parameter used to measure the proportion of high-fitness samples in the overall sample. It determines what proportion of high-fitness samples the algorithm prioritizes to guide the search direction or update the probability distribution in each iteration or update. A higher elite ratio can accelerate convergence, allowing the algorithm to focus on the current optimal region more quickly, but may lead to reduced search diversity and a higher risk of getting trapped in local optima. Conversely, a lower elite ratio retains the influence of more suboptimal samples, enhancing the exploratory nature of the search and the diversity of distributions, thus preventing premature convergence. By reasonably setting and dynamically adjusting the elite ratio, the algorithm can achieve a balance between global exploration and local refinement, improving optimization efficiency and inversion accuracy.

[0109] It should be noted that if the covariance of the first sampling distribution of cross-entropy is less than or equal to the covariance threshold, it indicates that the current sampling distribution is too concentrated, with most samples distributed in local regions, lacking sufficient exploration of the space. This state results in insufficient coverage of candidate samples, and the algorithm cannot effectively reach potential superior regions, thus reducing search diversity. Insufficient search diversity means that the optimization process relies too much on a few local superior samples. Particles or samples are easily concentrated near local optima and iterate repeatedly, lacking the ability to escape local regions, which may eventually lead to premature convergence.

[0110] The covariance deviation value is obtained by subtracting the covariance threshold from the first sampling distribution covariance of the cross-entropy. The covariance deviation value is divided by the first sampling distribution covariance of the cross-entropy to obtain the covariance deviation ratio. The covariance deviation ratio is multiplied by the preset maximum allowable reduction of the elite ratio to obtain the reduction of the elite ratio. The reduction of the elite ratio is added to the current elite ratio to obtain the updated elite ratio value.

[0111] Furthermore, the expansion covariance matrix is ​​analyzed using the following specific methods:

[0112] Calculate the mean square error of the samples within the sampling distribution obtained from the previous iteration, and then perform a weighted average of the samples based on their fitness to obtain a new expected vector.

[0113] Subtract the expected vector from each sample vector to obtain the bias vector, perform an outer product operation, and then sum the results according to their weights to form a new weighted variance estimate.

[0114] The original covariance matrix is ​​replaced with this weighted variance estimate or it is weighted and fused with it according to a certain learning rate, thereby realizing the covariance update.

[0115] ,

[0116] Where Σk+1 represents the covariance matrix of the next round, Σk represents the covariance matrix of the sampling distribution of the current round, β represents the learning rate, and Cov(Xe) represents the covariance matrix of the elite samples selected in the current round.

[0117] When the optimization phase is a local refinement phase, a multi-algorithm collaborative adjustment strategy is determined and implemented.

[0118] like Figure 4 The multi-algorithm collaborative adjustment strategy diagram for seismic parameter inversion based on the RBF surrogate model and multi-algorithm collaborative optimization includes the following steps: Continuously monitor the poor solution acceptance rate and the distribution covariance of cross-entropy sampling of the annealing algorithm, denoted as the second poor solution acceptance rate and the second sampling distribution covariance of the annealing algorithm, respectively. If the current optimal fitness has reached or exceeded the threshold, it is determined that no adjustment is needed, and the multi-algorithm collaboratively maintains the original state. If the optimal fitness has not reached the threshold, a dual judgment criterion is further activated. Based on the combination of these two criteria, a differentiated adjustment strategy is executed: if neither criterion is satisfied, the unadjusted state is maintained and iteration continues; if only the first judgment criterion is satisfied, the annealing temperature adjustment strategy is triggered; if only the second judgment criterion is satisfied, the cross-entropy distribution adjustment strategy is activated; when both criteria are satisfied simultaneously, it is determined that the system is in a local deadlock state, and the corresponding local deadlock integrated adjustment mechanism is immediately activated.

[0119] Furthermore, the multi-algorithm collaborative adjustment strategy is analyzed in the following ways:

[0120] In the later local refinement stage, the acceptance rate of inferior solutions of the annealing algorithm and the covariance of the cross-entropy sampling distribution are continuously monitored and denoted as the second inferior solution acceptance rate of the annealing algorithm and the second cross-entropy sampling distribution covariance, respectively.

[0121] Extract the preset fitness threshold from the database.

[0122] If the optimal fitness is greater than or equal to the fitness threshold, multi-algorithm collaborative adjustment is recorded as no adjustment.

[0123] It should be noted that if the optimal fitness is greater than or equal to the fitness threshold, it means that a sufficiently good solution has been found in the current search process, the algorithm is in the natural convergence stage, the optimal value of the objective function is close to or reaches the preset expectation, and it can continue to iterate without additional adjustment of the search strategy or parameters. The algorithm can maintain the current search state and continue to iterate in the good region, and stably converge to the optimal solution. Multi-algorithm collaborative adjustment is recorded as no adjustment.

[0124] If the optimal fitness is less than the fitness threshold, the second inferior solution acceptance rate of the annealing algorithm being lower than the inferior solution acceptance rate threshold is recorded as the first criterion for local refinement, and the second sampling distribution covariance of the cross-entropy being less than or equal to the covariance threshold is recorded as the second criterion for local refinement.

[0125] It should be noted that if the optimal fitness is less than the fitness threshold, it means that the optimal solution found in the current iteration has not yet reached the pre-set expected fitness level. This indicates that the search results are still some distance from the global optimum or the target requirement, and have not yet fully converged in the current state. There may still be better solutions in the search space that have not been explored. At the same time, the current samples or particles may be overly concentrated in local regions and lack coverage of potentially excellent global regions. The algorithm still has room for improvement and optimization. It is necessary to guide the samples or particles to explore the search space more effectively, improve search efficiency, avoid getting trapped in local optima, and thus gradually approach the global optimum.

[0126] If neither the first criterion for local refinement nor the second criterion for local refinement exists, the multi-algorithm collaborative adjustment is recorded as no adjustment, and the iterative search continues to generate the inversion execution signal.

[0127] It should be noted that if neither the first nor the second criterion for local refinement exists, it means that the acceptance rate of inferior solutions by the annealing algorithm is within a reasonable range in the current local refinement stage, and a moderate probability of accepting inferior solutions can be maintained. Furthermore, if the covariance of the cross-entropy sampling distribution is still greater than the threshold, it indicates that the sample distribution is not excessively concentrated and the search diversity is acceptable. There is no need to adjust the annealing temperature or the cross-entropy sampling distribution. Iterative optimization can continue under the current strategy to generate execution signals for inversion.

[0128] If the first criterion for local refinement exists and the second criterion for local refinement does not exist, the multi-algorithm collaborative adjustment is recorded as the execution annealing temperature adjustment strategy.

[0129] It should be noted that if the first criterion for local refinement exists but the second criterion does not, it means that during the local refinement stage, the acceptance rate of inferior solutions by the annealing algorithm is lower than the threshold, and the ability to accept inferior solutions needs to be improved to maintain search activity. However, if the covariance of the cross-entropy sampling distribution is still greater than the threshold, it indicates that the sampling distribution is not too concentrated and the search diversity remains normal. Therefore, the multi-algorithm collaborative adjustment strategy will implement the annealing temperature adjustment strategy. By appropriately adjusting the annealing temperature, the probability of the algorithm accepting inferior solutions will be increased, thereby activating local search and preventing it from getting stuck in local stagnation. At the same time, the distribution structure of the cross-entropy sampling will remain unchanged, providing more sufficient exploration capabilities for subsequent iterations.

[0130] If the first criterion for local refinement does not exist but the second criterion for local refinement does exist, the multi-algorithm collaborative adjustment is denoted as the implementation of the cross-entropy distribution adjustment strategy.

[0131] It should be noted that if the first criterion for local refinement does not exist and the second criterion for local refinement exists, it means that the acceptance rate of inferior solutions by the annealing algorithm is within a reasonable range, and there is no need to further increase the probability of accepting inferior solutions. However, if the cross-entropy sampling distribution covariance is less than or equal to the threshold, it indicates that the current sampling distribution is too concentrated and the search diversity is insufficient, which can easily lead to premature local convergence. Therefore, the multi-algorithm collaborative adjustment strategy will implement the cross-entropy distribution adjustment strategy, which expands the sampling range by expanding the covariance matrix and appropriately reduces the proportion of elites according to the covariance deviation value, thereby restoring search diversity, enhancing the coverage of unexplored areas around the local refinement region, and preventing local stagnation.

[0132] If both the first and second local refinement criteria exist, the multi-algorithm collaborative adjustment is denoted as local deadlock adjustment.

[0133] It should be noted that if both the first and second local refinement criteria exist, it means that during the local refinement stage, the acceptance rate of inferior solutions by the annealing algorithm is lower than the threshold, and the covariance of the cross-entropy sampling distribution is also less than or equal to the threshold. This indicates that the algorithm currently lacks the ability to accept inferior solutions and the sampling distribution is too concentrated, resulting in a serious lack of search diversity. The algorithm may fall into a local deadlock state and find it difficult to continue to effectively explore the optimal solution region. In this case, the multi-algorithm collaborative adjustment strategy will implement local deadlock adjustment.

[0134] Further, local lockout adjustment, the specific analysis method is as follows:

[0135] The annealing temperature increase and freezing round ratio are determined by subtracting the second inferior solution acceptance rate threshold from the inferior solution acceptance rate of the annealing algorithm. The annealing temperature freezing round is then determined based on the freezing round ratio and the preset maximum freezing round.

[0136] It should be noted that in this embodiment, a mapping table or function model is pre-established between the inferior solution acceptance rate deviation value, the annealing temperature increase value, and the freezing round ratio, and stored in the configuration file for unified algorithm invocation. The preset logic of the mapping table is as follows: when the inferior solution acceptance rate deviation value is large, it corresponds to a higher annealing temperature increase value and a longer freezing round ratio, so as to significantly increase the probability of accepting inferior solutions and ensure temperature stability in multiple iterations after temperature increase, thereby activating search activity and expanding the exploration range; when the deviation value is small or close to the target value, it corresponds to a smaller temperature increase value and a shorter freezing round ratio, so as to avoid excessive perturbation and ensure the stability of local refinement search. In actual operation, the optimization algorithm first calculates the second inferior solution acceptance rate of the current iterative annealing algorithm minus the inferior solution acceptance rate threshold to obtain the inferior solution acceptance rate deviation value. Then, the inferior solution acceptance rate deviation value is input into the mapping table or function model to find or calculate the corresponding annealing temperature increase value and freezing round ratio, and applied to the annealing temperature update and freezing strategy, thereby realizing the adaptive control of temperature increase and freezing round, enabling the algorithm to dynamically restore search activity and promote optimization iteration in a local dead state.

[0137] The freezing round ratio is multiplied by the preset maximum freezing round to obtain the annealing temperature freezing round.

[0138] The annealing temperature update value is obtained by adding the annealing temperature increase value to the current annealing temperature recorded in the system log. The freezing round ratio is multiplied by the preset maximum freezing round to obtain the annealing temperature freezing round, which is used to keep the temperature stable in the next few iterations, so that the algorithm maintains a high search activity within a fixed period and prevents the temperature from dropping too quickly until the freezing round ends. The annealing temperature remains unchanged within the freezing round.

[0139] The number of random samples is determined based on the deviation between the covariance of the second sampling distribution of cross-entropy and the covariance threshold. The proportion of new random samples injected into the neighborhood and global range of the current cross-entropy distribution is determined based on the deviation between the optimal fitness and the fitness threshold.

[0140] It should be noted that the covariance bias value is obtained by subtracting the covariance threshold from the second sampling distribution covariance of the cross-entropy, and the number of random samples is determined accordingly.

[0141] In this embodiment, a mapping table or function model between the covariance deviation value and the number of randomly generated samples is pre-established and stored in a configuration file for unified algorithm invocation. The preset logic of this mapping table is as follows: when the covariance deviation value is large, it indicates that the current sampling distribution is relatively dispersed and the search diversity is sufficient, and the number of randomly generated samples can be appropriately reduced to reduce computational overhead; when the covariance deviation value is small or close to the threshold, it indicates that the sampling distribution is too concentrated and the search diversity is insufficient, and more random samples should be generated to expand the sample coverage, improve the exploration capability, and prevent local convergence. In actual operation, the optimization algorithm first calculates the covariance of the cross-entropy sampling distribution of the current iteration and compares it with the covariance threshold to obtain the covariance deviation value; then, the deviation value is input into the mapping table or function model to find or calculate the corresponding number of random samples, and the corresponding number of random samples are generated in the next iteration, thereby realizing adaptive control of sample generation, enabling the algorithm to dynamically adjust the exploration intensity and coverage under different search states.

[0142] The fitness bias is obtained by subtracting the fitness threshold from the optimal fitness. This bias determines the proportion of new random samples injected into the current cross-entropy distribution neighborhood and the global scope. In this embodiment, a mapping table or function model between the fitness bias and the proportion of new random samples injected is pre-established and stored in a configuration file for unified algorithm calls. The preset logic of this mapping table is as follows: when the fitness bias is large, a higher proportion of new samples is injected to introduce more random samples into the current sampling distribution neighborhood and the global scope, thereby enhancing search diversity, preventing local deadlock, and improving the exploration capability in the local refinement stage; when the bias is small or close to the threshold, a lower or zero proportion of new samples is injected to maintain the existing sample distribution structure and search stability. In actual operation, the optimization algorithm first calculates the fitness bias between the optimal fitness and the fitness threshold of the current iteration; then, the fitness bias is input into the mapping table or function model to find or calculate the corresponding proportion of new samples injected, and random samples are generated in the cross-entropy distribution neighborhood and the global scope according to this proportion for the next iteration, so as to dynamically adjust the sample coverage and search activity, and realize adaptive diversity control in the local refinement stage.

[0143] At the same time, a short-term taboo mechanism is activated to extract regions in the system log that continuously exhibit local deadlock adjustments. The sampling weight reduction value is determined based on the number of local deadlock adjustments in that region, guiding the search to migrate to unexplored regions.

[0144] It should be noted that regions with consecutive local deadlock adjustments are extracted from the system logs, and the number of times each region triggers a local deadlock adjustment is counted. In this embodiment, the number of local deadlocks in consecutively failing regions is first counted. Then, this number is used as input to search or calculate by referring to a pre-established weight reduction mapping table or applying preset weight reduction rules. Each deadlock frequency interval in the mapping table corresponds to a sampling weight reduction ratio, or the rule maps the frequency to the weight reduction magnitude through linear, non-linear, or exponential functions. In this way, the algorithm dynamically determines the sampling weight that should be reduced for the region in the next iteration based on the number of deadlocks. The higher the frequency, the greater the corresponding weight reduction value, thereby ensuring that the sampling attention of consecutively failing regions is gradually reduced, while freeing up more sampling opportunities for unexplored regions, maintaining overall search diversity and optimization efficiency.

[0145] After receiving the inversion execution signal, a closed-loop optimization inversion is performed to obtain the globally optimal inversion result.

[0146] Seismic parameter inversion methods combine an RBF surrogate model with particle swarm optimization, annealing, and cross-entropy sampling algorithms. During the iterative process, candidate formation parameter samples are generated, and the surrogate model rapidly predicts their corresponding observed responses. Then, forward modeling is used to realistically simulate key samples, and the simulated waveforms are compared with observed waveforms to guide the algorithm in updating its search direction. In this way, with the assistance of the surrogate model, the three algorithms combine global search with local refinement, gradually approaching the optimal formation parameter solution, thus completing the joint inversion of P-wave velocity, S-wave velocity, and density, and obtaining the globally optimal inversion result.

[0147] The above description is only an optional embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

Claims

1. A seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization, characterized in that, The method includes: Raw seismic observation data is collected, standardized, and then an RBF surrogate model is constructed for inversion and iterative optimization. During the inversion iteration process, the timestamps of the shared states between algorithms are collected in real time, proxy verification is performed, the adjustment requirements of the RBF proxy model are determined, and the optimization stage is dynamically divided. When the optimization phase is a global exploration phase, determine the annealing temperature adjustment strategy and the cross-entropy distribution adjustment strategy. When the optimization phase is a local refinement phase, a multi-algorithm collaborative adjustment strategy is determined and implemented. After receiving the inversion execution signal, a closed-loop optimization inversion is performed to obtain the globally optimal inversion result.

2. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 1, characterized in that, The specific analysis method for constructing the RBF proxy model is as follows: The raw seismic observation data is standardized, and the reasonable range of inversion parameters is determined based on prior geological knowledge. The Latin hypercube sampling technique is used to generate representative initial sample points that uniformly cover the entire parameter space. Forward modeling is performed on each initial sample point to evaluate its fit with the observed data; The best-fitting sample point is selected from the initial sample points as the center point of the radial basis function model. The RBF surrogate model is constructed using cubic spline basis functions. The optimal weight parameters of the RBF surrogate model are determined by solving the linear equation system.

3. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 1, characterized in that, The execution proxy verification determines the adjustment requirements of the RBF proxy model. The specific analysis method is as follows: Each initial sample point is used to construct an initial sample set, which is then used for iterative optimization using multiple algorithms. The timestamps of shared states between algorithms are collected in real time. When it is detected that the timestamp read by any algorithm is inconsistent with the current timestamp, the cross-algorithm data interaction of the corresponding algorithm is temporarily suspended, and proxy verification is performed first. The specific analysis process is as follows: The current iteration stratigraphic parameter vector is obtained as a sample and input into the forward modeling module for simulation to obtain the simulated seismic waveform. The simulated seismic waveform is compared with the observed seismic waveform, and the mean square error is calculated. Extract the preset mean square error threshold from the database; If the mean square error is greater than or equal to the mean square error threshold, the adjustment requirement of the RBF proxy model is recorded as a demand adjustment. The RBF kernel width is adjusted to correct the proxy model by the deviation between the mean square error and the mean square error threshold, and the local state cache is refreshed. If the mean squared error is less than the mean squared error threshold, the adjustment requirement of the RBF proxy model is marked as no adjustment is needed, the use of this sample for cross-algorithm updates is temporarily suspended, and a local state refresh is triggered.

4. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 1, characterized in that, The specific analysis method for the dynamic partitioning and optimization stage is as follows: Calculate the rate of change of the optimal fitness in each iteration: Extract the preset optimal fitness change rate threshold from the database; When the rate of change of the optimal fitness is greater than or equal to the rate of change threshold, the optimization phase is divided into the early global exploration phase. When the rate of change of the optimal fitness is less than the rate of change threshold, the optimization phase is divided into the late local refinement phase.

5. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 1 further includes adaptive adjustment of PSO algorithm parameters, the specific analysis method of which is as follows: Extract the preset initial inertia weight, final inertia weight, and maximum number of iterations from the database; Based on the ratio of the current iteration count to the maximum iteration count, the current inertia weight is determined relative to the initial and final inertia weights, thus determining the current inertia weight of the PSO particle swarm.

6. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 1, characterized in that, The specific analysis methods for determining the annealing temperature adjustment strategy and the cross-entropy distribution adjustment strategy are as follows: In the early global exploration phase, the acceptance rate of the inferior solution of the annealing algorithm is denoted as the first inferior solution acceptance rate of the annealing algorithm, and the covariance of the cross-entropy sampling distribution is denoted as the first cross-entropy sampling distribution covariance. Extract the pre-defined poor solution acceptance rate threshold and covariance threshold from the database; If the acceptance rate of the first inferior solution in the annealing algorithm is lower than the inferior solution acceptance rate threshold, it is recorded as the first condition for global adjustment. If the covariance of the first sampling distribution of cross-entropy is less than or equal to the covariance threshold, it is denoted as the second condition for global adjustment. If there is a first condition for global adjustment and there is no second condition for global adjustment, then the global exploration phase adjustment strategy is recorded as the execution annealing temperature adjustment strategy. If there is a second global adjustment condition but no first global adjustment condition, then the global exploration phase adjustment strategy is recorded as executing the cross-entropy distribution adjustment strategy. If neither of the two conditions exists, the global exploration phase adjustment strategy is recorded as continuing iterative optimization, and an inversion execution signal is generated.

7. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 6, characterized in that, The specific analysis method for implementing the annealing temperature adjustment strategy is as follows: The deviation value of the acceptance rate of the first inferior solution in the annealing algorithm is obtained based on the acceptance rate of the first inferior solution and the threshold value of the acceptance rate of the inferior solution. The reduction value of the annealing temperature decrease rate is determined based on the deviation value of the acceptance rate of the inferior solution. The acceptance rate of the inferior solution in the annealing algorithm is statistically analyzed in real time within the sliding window. The temperature decay coefficient is corrected based on the reduction value of the annealing temperature decrease rate. When the acceptance rate of inferior solutions in the adjusted annealing algorithm is still lower than the threshold for the acceptance rate of inferior solutions, the algorithm is rolled back to the previous stable temperature and the superimposed short-term disturbance factor is determined based on the quadratic deviation between the acceptance rate of inferior solutions in the adjusted annealing algorithm and the threshold for the acceptance rate of inferior solutions, so as to restore search activity. The specific analysis method for the cross-entropy distribution adjustment strategy is as follows: The expansion of the covariance matrix and the reduction of the elite ratio are determined by combining the deviation of the first sampling distribution covariance of the cross-entropy and the covariance threshold with the preset maximum allowable reduction of the elite ratio, thereby obtaining the updated elite ratio value, which is used in the next iteration.

8. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 7, characterized in that, The specific analysis method for the dilated covariance matrix is ​​as follows: Calculate the objective function value of the samples within the sampling distribution obtained from the previous iteration, and weight the samples according to their fitness to obtain a new expectation vector; Perform an outer product operation on the deviation vector between each sample and the expected vector, and then sum them according to their weights to form a new weighted variance estimate; The original covariance matrix is ​​replaced with this weighted variance estimate or it is weighted and fused with it according to a certain learning rate, thereby realizing the covariance update.

9. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 1, characterized in that, The specific analysis method for the multi-algorithm collaborative adjustment strategy is as follows: In the later local refinement stage, the acceptance rate of inferior solutions of the annealing algorithm and the covariance of the cross-entropy sampling distribution are continuously monitored and denoted as the second inferior solution acceptance rate of the annealing algorithm and the second cross-entropy sampling distribution covariance, respectively. Extract the preset fitness threshold from the database; If the optimal fitness is greater than or equal to the fitness threshold, multi-algorithm collaborative adjustment is recorded as no adjustment; If the optimal fitness is less than the fitness threshold, the second inferior solution acceptance rate of the annealing algorithm being lower than the inferior solution acceptance rate threshold is recorded as the first criterion for local refinement, and the second sampling distribution covariance of the cross-entropy being less than or equal to the covariance threshold is recorded as the second criterion for local refinement. If neither the first criterion for local refinement nor the second criterion for local refinement exists, the multi-algorithm collaborative adjustment is recorded as no adjustment, and the iterative search continues to generate the inversion execution signal; If the first criterion for local refinement exists and the second criterion for local refinement does not exist, the multi-algorithm collaborative adjustment is recorded as the execution annealing temperature adjustment strategy. If the first criterion for local refinement does not exist but the second criterion for local refinement exists, the multi-algorithm collaborative adjustment is denoted as the implementation of the cross-entropy distribution adjustment strategy. If both the first and second local refinement criteria exist, the multi-algorithm collaborative adjustment is denoted as local deadlock adjustment.

10. The seismic parameter inversion method based on the RBF surrogate model and multi-algorithm collaborative optimization as described in claim 9, characterized in that, The specific analysis method for the local lockout adjustment is as follows: The annealing temperature increase and freezing round ratio are determined based on the deviation between the second inferior solution acceptance rate and the inferior solution acceptance rate threshold of the annealing algorithm. The annealing temperature freezing round is determined based on the freezing round ratio and the preset maximum freezing round. The number of random samples is determined based on the deviation between the covariance of the second sampling distribution of cross-entropy and the covariance threshold. The proportion of new random samples injected into the neighborhood and global range of the current cross-entropy distribution is determined based on the deviation between the optimal fitness and the fitness threshold. At the same time, a short-term taboo mechanism is activated to extract regions in the system log that continuously exhibit local deadlock adjustments. The sampling weight reduction value is determined based on the number of local deadlock adjustments in that region, guiding the search to migrate to unexplored regions.

Images