Spherical wave oscillation-slowly-varying kernel high-order decoupling method and device for physical property inversion

Abstract

The application relates to a spherical wave oscillation-slowly-varying kernel high-order decoupling method and device facing physical property inversion, and belongs to the field of geophysical exploration and acoustic exploration. The method comprises the following steps: obtaining observation data and inputting, initializing seabed elastic parameters; determining an integral interval according to the observation data, decomposing the integral interval into a plurality of subintervals and setting a decoupling order; calculating local integral matrixes of each order of oscillation kernel corresponding to the incidence angle of the observation data; calculating angle partial derivative vectors of each order of slowly-varying kernel; calculating each order of spherical wave reflection coefficient vectors; calculating high-order decoupled spherical wave reflection coefficients; and solving a subsequent inversion problem by using the obtained high-order decoupled spherical wave reflection coefficients. The application effectively improves the calculation efficiency of seabed spherical wave reflection coefficients in the inversion process; improves the calculation precision by increasing the decoupling order, or reduces the number of subintervals by increasing the decoupling order under the condition of certain calculation precision.

Description

Technical Field

[0001] This invention belongs to the field of geophysics, specifically relating to a high-order decoupling method and apparatus for spherical wave oscillation-slowly varying kernels for property inversion. Background Technology

[0002] The relationship between seabed reflection coefficient and incident angle carries rich information about the physical properties of seabed sediments. Due to its simple and efficient calculation process, plane wave theory has been widely used in seabed reflection effect analysis and physical property inversion. However, in actual geophysical exploration or acoustic detection, seismic waves are spherical waves generated by point sources. Plane waves are an approximation of spherical waves in the high-frequency, far-field case. This means that plane wave theory cannot accurately describe the seabed spherical wave reflection coefficient in the low-frequency, near-field case, especially at large incident angles where the plane wave reflection coefficient will produce significant errors. While spherical wave theory can better describe the relationship between seabed reflection coefficient and incident angle, its computational efficiency is significantly lower than that of calculation methods based on plane wave theory because it requires integration of rapidly oscillating functions. This severely limits the application of spherical wave theory in related forward and inverse problems.

[0003] The seabed spherical wave reflection coefficient can be expressed as a weighted superposition of plane wave reflection coefficients at different incident angles. The weighting coefficients are determined by the wave field frequency, observation system parameters, and seawater properties, and oscillate rapidly with the incident angle. In contrast, the plane wave reflection coefficient is determined only by the properties of the seabed on both sides, and its change with the incident angle is slower compared to the weighting coefficients. The rapidly oscillating weighting coefficients constitute the oscillation kernel of the spherical wave reflection coefficient, while the slowly changing plane wave reflection coefficient constitutes the gradually varying kernel. Based on the aforementioned characteristics of the oscillation kernel and the gradually varying kernel, this invention performs a Taylor expansion of the gradually varying kernel within local sub-intervals, thereby achieving arbitrary-order decoupling of the oscillation and gradually varying kernels. Arbitrary-order decoupling achieves a balance between sub-interval discretization and computational accuracy. With a fixed number of local sub-intervals, computational accuracy can be improved by increasing the decoupling order; conversely, with a fixed computational accuracy, the number of sub-intervals can be reduced by increasing the decoupling order. The decoupling of the oscillation-gradient kernel allows for the calculation and updating of only the gradient kernel during the inversion process, without the need to repeatedly calculate the oscillation kernel. This avoids the repeated integration calculation of the oscillation function in the traditional inversion process and significantly improves the calculation efficiency of the spherical wave reflection coefficient during the inversion process. Summary of the Invention

[0004] By decoupling the oscillating kernel and the slowly varying kernel in the integral expression of the spherical wave reflection coefficient, the redundant calculation of the oscillating kernel part during the inversion process can be avoided, effectively improving the calculation efficiency of the seabed spherical wave reflection coefficient during the inversion process. Furthermore, higher-order decoupling methods can achieve a balance between sub-interval discretization and computational accuracy by varying the decoupling order; with a fixed number of local sub-intervals, increasing the decoupling order can improve computational accuracy, or with a fixed computational accuracy, increasing the decoupling order can reduce the number of sub-intervals.

[0005] This invention is achieved through the following technical solution:

[0006] A method for decoupling spherical wave oscillations and slowly varying kernels for property inversion, comprising the following steps:

[0007] Step 1: Obtain observation data And input, initialize the seabed elastic parameters ;

[0008] Step 2: Based on the observation data Determine the integration interval, decompose the integration interval into several sub-intervals, and set the decoupling order;

[0009] Furthermore, in this step, based on the observation data input in the first step... Determine the integration interval Undetermined parameters in The integration interval is discretized into several sub-intervals, and then... Represents any one of the subintervals ( ). The value of does not need to be strictly calculated; it is only necessary to ensure that the integral outside the interval is much smaller than the integral inside the interval. The union of all subintervals forms the integration interval. The intersection of any two subintervals is an empty set. Let the decoupling order be denoted as . .

[0010] Step 3: Calculate the local integral matrix of each order of oscillation kernel corresponding to the incident angle of the observed data;

[0011] In this step, the local integral matrix of each order of oscillation kernel is used It means that among them , representing zeroth order, first order, second order, ... Order; Local integral matrix of each order of oscillation kernel ;in, ... Indicates the angle of incidence corresponding to the observed data. ; Indicates a subinterval under a given incident angle. Internal oscillating nucleus of The specific formula for calculating the value of the first-order definite integral is as follows:

[0012] ,in The wave number of the seismic wave. Indicates the propagation distance of the wave field. It is a natural constant. express factorial, sub-interval Taylor series expansion points within;

[0013] Step 4: Calculate the angular deflection parameters of each order of the slowly varying kernel;

[0014] In sequence, in all sub-intervals of "Step Two" Internally determined by the current seabed elastic parameters Calculate the plane wave reflection coefficient For the angle of incidence The partial derivatives of each order yield the angular partial derivatives of the slowly varying kernel. ;in, , , The meaning is the same as in "Step Three".

[0015] When calculating the angular deflection quantities of each order of the slowly varying core for the first time, the current seabed elastic parameters... Refers to the seabed elastic parameters initialized in the "first step". In subsequent calculations of the angular deflection derivatives of the slowly varying core, the current seabed elastic parameters are determined during the iterative solution process. Refers to the updated seabed elasticity parameters after "step seven". .

[0016] Step 5: Calculate the reflection coefficient vectors of each order of spherical wave; in this step, the first... First-order reflection coefficient vector ,in, and It is obtained from "step three" and "step four" respectively.

[0017] Step 6: Calculate the reflection coefficient of higher-order decoupled spherical waves. .

[0018] Step 7: Use the higher-order decoupled spherical wave reflection coefficients obtained in Step 6 to solve the subsequent inversion problem. For example... Figure 1 As shown in the shaded area, in the subsequent inversion problem solution process, it is only necessary to repeat "steps four, five and six" after updating the seabed elastic parameters, without having to repeat the calculation of the local integral matrices of each order of oscillation kernel in step three.

[0019] This invention also provides a high-order decoupling system for seismic spherical reflection wave oscillation-gradient kernel for seabed physical parameter inversion; the system includes a data input and initialization module, an integral interval decomposition module, a module for calculating the local integral matrix of each order oscillation kernel, a module for calculating the angular deflection vectors of each order of the gradually varying kernel, a module for calculating the spherical wave reflection coefficient vectors of each order, a module for calculating the spherical wave reflection coefficients of the high-order decoupled system, and a solution module for the subsequent inversion problem;

[0020] The data input module and initialization module execute the first step of the method; the integral interval decomposition module executes the second step; the module for calculating the local integral matrix of each oscillating kernel executes the third step; the module for calculating the angular deflection vectors of each order of the slowly varying kernel executes the fourth step; the module for calculating the spherical wave reflection coefficient vectors of each order executes the fifth step; the module for calculating the higher-order decoupled spherical wave reflection coefficient executes the sixth step; and the subsequent inversion problem solution module executes the sixth step. The specific steps are as follows:

[0021] The data input module and initialization module are used to input observation data. Initialize seabed elastic parameters ;

[0022] The integral interval decomposition module is used to decompose the integral interval into several sub-intervals and set the decoupling order;

[0023] The module for calculating the local integral matrix of each oscillation kernel is used to calculate the local integral matrix of each oscillation kernel corresponding to the incident angle of the observed data.

[0024] The module for calculating the angular deflection parameters of each order of the slowly varying kernel is used to calculate the angular deflection parameters of each order of the slowly varying kernel; it sequentially decomposes all sub-intervals of the integration interval module. Internally determined by the current seabed elastic parameters Calculate the plane wave reflection coefficient For the angle of incidence The partial derivatives of each order yield the angular partial derivatives of the slowly varying kernel. ;in, , ;

[0025] The module for calculating the reflection coefficient vector of each order of spherical wave is used to calculate the reflection coefficient vector of each order of spherical wave; in this step, the first... First-order reflection coefficient vector ,in and The results are obtained by the module for calculating the local integral matrix of each order of oscillating kernel and the module for calculating the angular deflection vectors of each order of slowly varying kernel, respectively.

[0026] The module for calculating the higher-order decoupled spherical wave reflection coefficient is used to calculate the higher-order decoupled spherical wave reflection coefficient. ;

[0027] The subsequent inversion problem solving module is used to solve the subsequent inversion problem by utilizing the higher-order decoupled spherical wave reflection coefficients in the higher-order decoupled spherical wave reflection coefficient calculation module. In the process of solving the subsequent inversion problem, it is only necessary to repeatedly execute the modules for calculating the angle deflection vectors of the slowly varying kernel, the modules for calculating the spherical wave reflection coefficient vectors of each order, and the modules for calculating the spherical wave reflection coefficient vectors of each order after updating the seabed elastic parameters.

[0028] The present invention also provides a high-order decoupling device for seismic spherical reflection wave oscillation-gradiently varying kernel for seabed physical parameter inversion, wherein the device is equipped with the above-mentioned decoupling system and performs the above-mentioned decoupling method.

[0029] The beneficial effects of this invention compared to existing technologies are as follows: By decoupling the oscillating kernel and the slowly varying kernel in the integral expression of the spherical wave reflection coefficient, the repeated calculation of the oscillating kernel part during the inversion process can be avoided, effectively improving the calculation efficiency of the seabed spherical wave reflection coefficient during the inversion process. Figure 1 (Shaded area). Furthermore, higher-order decoupling methods can achieve a balance between sub-interval discretization and computational accuracy by varying the decoupling order; computational accuracy can be improved by increasing the decoupling order when the number of local sub-intervals is fixed, or the number of sub-intervals can be reduced by increasing the decoupling order when computational accuracy is fixed. Attached Figure Description

[0030] Figure 1 Flowchart of a high-order decoupling method for seismic spherical reflection wave oscillation-gradient kernel for seabed physical parameter inversion;

[0031] Figure 2 The following is a diagram of observation data for an example: A represents low-velocity seabed data, and B represents high-velocity seabed data.

[0032] Figure 3 This is a diagram of the local integral matrix of the oscillating kernel. A is the natural logarithm of the absolute value of the real part of the zeroth order, B is the natural logarithm of the absolute value of the imaginary part of the zeroth order, C is the natural logarithm of the absolute value of the real part of the first order, D is the natural logarithm of the absolute value of the imaginary part of the first order, E is the natural logarithm of the absolute value of the real part of the second order, and F is the natural logarithm of the absolute value of the imaginary part of the second order.

[0033] Figure 4 This is a diagram of the local integral matrix of the oscillating kernel. A is the natural logarithm of the absolute value of the real part of the zeroth order, B is the natural logarithm of the absolute value of the imaginary part of the zeroth order, C is the natural logarithm of the absolute value of the real part of the first order, D is the natural logarithm of the absolute value of the imaginary part of the first order, E is the natural logarithm of the absolute value of the real part of the second order, and F is the natural logarithm of the absolute value of the imaginary part of the second order.

[0034] Figure 5 This is a directional diagram of the angle deflection of the slowly varying core at low speeds on the seabed. A is the zeroth-order real part, B is the zeroth-order imaginary part, C is the first-order real part, D is the first-order imaginary part, E is the second-order real part, and F is the second-order imaginary part.

[0035] Figure 6 This is a directional diagram of the angle deflection of the slowly varying core at low speeds on the seabed. A is the zeroth-order real part, B is the zeroth-order imaginary part, C is the first-order real part, D is the first-order imaginary part, E is the second-order real part, and F is the second-order imaginary part.

[0036] Figure 7 This is a directional vector diagram showing the angle deflection of different orders of hypersonic, slowly varying seabed nuclei. A is the zeroth-order real part, B is the zeroth-order imaginary part, C is the first-order real part, D is the first-order imaginary part, E is the second-order real part, and F is the second-order imaginary part.

[0037] Figure 8 This is a directional vector diagram showing the angle deflection of different orders of hypersonic, slowly varying seabed nuclei. A is the zeroth-order real part, B is the zeroth-order imaginary part, C is the first-order real part, D is the first-order imaginary part, E is the second-order real part, and F is the second-order imaginary part.

[0038] Figure 9 This is a comparison chart showing the relationship between the seabed reflection coefficient at low speeds of sound and the incident angle. In A, the black solid line represents the zero-order decoupled reflection coefficient obtained by the present invention; the green solid line represents the second-order decoupled reflection coefficient obtained by the present invention; and the red dot in A represents the accurate reflection coefficient obtained by equation (1). In B, the black solid line is the difference between the zero reflection coefficient and the accurate reflection coefficient of the present invention; and the green solid line is the difference between the second-order reflection coefficient and the accurate reflection coefficient of the present invention.

[0039] Figure 10 This is a comparison chart showing the relationship between the seabed reflection coefficient at low speeds of sound and the incident angle. In A, the black solid line represents the zero-order decoupled reflection coefficient obtained by the present invention; the green solid line represents the second-order decoupled reflection coefficient obtained by the present invention; and the red dot in A represents the accurate reflection coefficient obtained by equation (1). In B, the black solid line is the difference between the zero reflection coefficient and the accurate reflection coefficient of the present invention; and the green solid line is the difference between the second-order reflection coefficient and the accurate reflection coefficient of the present invention.

[0040] Figure 11 This is a comparison chart showing the relationship between the hypersonic seabed reflection coefficient and the incident angle. In A, the black solid line represents the zero-order decoupled reflection coefficient obtained by the present invention; the green solid line represents the second-order decoupled reflection coefficient obtained by the present invention; and the red dot in A represents the accurate reflection coefficient obtained by equation (1). In B, the black solid line is the difference between the zero reflection coefficient and the accurate reflection coefficient of the present invention; and the green solid line is the difference between the second-order reflection coefficient and the accurate reflection coefficient of the present invention.

[0041] Figure 12 This is a comparison chart showing the relationship between the hypersonic seabed reflection coefficient and the incident angle. In A, the black solid line represents the zero-order decoupled reflection coefficient obtained by the present invention; the green solid line represents the second-order decoupled reflection coefficient obtained by the present invention; and the red dot in A represents the accurate reflection coefficient obtained by equation (1). In B, the black solid line is the difference between the zero reflection coefficient and the accurate reflection coefficient of the present invention; and the green solid line is the difference between the second-order reflection coefficient and the accurate reflection coefficient of the present invention.

[0042] Figure 13 The hypersonic seabed posterior probability density map is generated by the exact equation-driven inversion. A is the joint probability distribution of density and P-wave velocity, B is the joint probability distribution of S-wave velocity and P-wave velocity, C is the joint probability distribution of S-wave velocity and density, D is the probability distribution of P-wave velocity, E is the probability distribution of density, and F is the probability distribution of S-wave velocity.

[0043] Figure 14 The diagram shows the posterior probability distribution of hypersonic seabed velocity after zero-order decoupling inversion, with I=450. A is the joint probability distribution of density and P-wave velocity, B is the joint probability distribution of S-wave velocity and P-wave velocity, C is the joint probability distribution of S-wave velocity and density, D is the probability distribution of P-wave velocity, E is the probability distribution of density, and F is the probability distribution of S-wave velocity.

[0044] Figure 15 The second-order decoupled inversion hypersonic seabed posterior probability density map is shown, I=450. A is the joint probability distribution of density-P-wave velocity, B is the joint probability distribution of S-wave velocity-P-wave velocity, C is the joint probability distribution of S-wave velocity-density, D is the probability distribution of P-wave velocity, E is the probability distribution of density, and F is the probability distribution of S-wave velocity.

[0045] Figure 16 The diagram shows the posterior probability distribution of hypersonic seabed velocity after zero-order decoupling inversion, where I=1500. A is the joint probability distribution of density and P-wave velocity, B is the joint probability distribution of S-wave velocity and P-wave velocity, C is the joint probability distribution of S-wave velocity and density, D is the probability distribution of P-wave velocity, E is the probability distribution of density, and F is the probability distribution of S-wave velocity.

[0046] Figure 17The diagram shows the posterior probability distribution of hypersonic seabed velocity in the second-order decoupled inversion, where I=1500. A is the joint probability distribution of density and P-wave velocity, B is the joint probability distribution of S-wave velocity and P-wave velocity, C is the joint probability distribution of S-wave velocity and density, D is the probability distribution of P-wave velocity, E is the probability distribution of density, and F is the probability distribution of S-wave velocity. Detailed Implementation

[0047] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings. It is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the protection scope of the present invention.

[0048] This invention relates to a high-order decoupling method for seismic spherical reflection wave oscillations and slowly varying kernels for seabed physical parameter inversion, such as... Figure 1 As shown, its main principle is as follows.

[0049] As shown in equation (1), the spherical wave reflection coefficient can be expressed as a weighted superposition of the plane wave reflection coefficients at different incident angles:

[0050] (1)

[0051] In equation (1), The plane wave reflection coefficient, These are the weighting coefficients. The wave number of the seismic wave. Indicates the propagation distance of the wave field. It is a natural constant. The plane wave reflection coefficient, The horizontal distance between the seismic source and the receiving point. and These represent the distances from the seismic source and the receiving point to the reflecting interface, respectively. Denotes the zeroth-order Bessel function of the first kind. Let the integration variable represent the integration interval. The incident angle inside. The accurate spherical wave reflection coefficient can be obtained from equation (1). However, the calculation process requires integration of a rapidly oscillating function, which is inefficient. Weighting coefficients The oscillating nucleus of the spherical wave reflection coefficient is determined by the wave field frequency, observation system parameters, and seawater properties, and oscillates rapidly with the incident angle. The plane wave reflection coefficient... Determined solely by the physical properties on both sides of the seabed, and by the weighting coefficients. In contrast, the oscillating kernel changes slowly with the incident angle, forming a slowly varying kernel for the spherical wave reflection coefficient. This invention decouples the oscillating kernel and the slowly varying kernel, avoiding redundant calculations of the oscillating kernel during the inversion process and effectively improving the computational efficiency of the seabed spherical wave reflection coefficient during both forward and inversion processes. Higher-order decoupling methods can achieve a balance between sub-interval discretization and computational accuracy by varying the decoupling order; computational accuracy can be improved by increasing the decoupling order when the number of local sub-intervals is fixed, or the number of sub-intervals can be reduced by increasing the decoupling order when computational accuracy is fixed.

[0052] The proposed method for seismic spherical reflection wave oscillation-gradiently varying kernel high-order decoupling for seabed physical parameter inversion is as follows: Figure 1 As shown, the specific steps include:

[0053] Step 1: Input observation data Initialize seabed elastic parameters .

[0054] This embodiment uses single-reflection point data under both low-speed and high-speed seabed conditions as examples to describe the technical solution. The observation data is as follows: Figure 2 As shown, A represents low-velocity seabed data, and B represents high-velocity seabed data; the minimum incident angle of the data is 1°, and the maximum incident angle is 60°, distributed at 1° intervals. The frequency of all data shown is 100Hz, and the distance from the source and receiver to the reflecting interface is 50m (i.e.,...). ), and It can be calculated from the incident angle based on geometric relationships. Elastic parameters of the seabed under two different seabed conditions. All initialized to .

[0055] Step 2: Divide the integration interval Decompose into several sub-intervals, using Represents any one of the subintervals ( The subinterval decomposition can be either equally spaced or non-equally spaced; the computational accuracy increases as the spacing decreases. The decoupling order is set. For the same spacing, the higher the decoupling order, the higher the calculation accuracy.

[0056] exist The decay is rapid at higher frequencies, and the decay is even faster at higher frequencies. The value of is chosen such that the integral outside the interval is much smaller than that inside the interval when the frequency is highest.

[0057] This embodiment takes The value is 1. To facilitate comparison of the impact of the number of subintervals on computational accuracy, the integration interval is decomposed into 450 and 1500 subintervals respectively, i.e. or Sub-interval index and Figures 3 to 8 The horizontal axis corresponds accordingly. To compare the impact of different decoupling orders on computational accuracy, this example uses the decoupling order... Set them to 0 and 2 respectively.

[0058] Step 3: Calculate the local integral matrix of each order of oscillation kernel corresponding to the incident angle of the observed data.

[0059] In this step, the local integral matrix of each order of oscillation kernel is used It means that among them , representing zeroth order, first order, second order, and so on. Local integral matrix of the oscillating kernel. .in, ... Indicates the angle of incidence corresponding to the observed data. . Indicates a subinterval under a given incident angle. Internal oscillating nucleus of The specific formula for calculating the value of the first-order definite integral is as follows:

[0060] .

[0061] The local integral matrix of the oscillation kernel calculated in this embodiment based on the input data and related parameters in "Step 1" and "Step 2" is as follows: Figure 3 , Figure 4 As shown. Figure 3 , Figure 4 The vertical axis represents the angle of incidence. ... and Figure 2 The input data shown corresponds to the incident angles, which are 1°, 2°, ..., 65° (n=65); the horizontal axis represents the sub-interval number, which corresponds to the sub-interval number in "Step Two". correspond.

[0062] Step 4: Calculate the angular deflection conductance of each order of the slowly varying kernel. This is done sequentially in all subintervals of "Step 2". Internal calculation of plane wave reflection coefficient The angular partial derivatives of each order are used to obtain the angular partial derivatives of the slowly varying kernel. .in, , , The meaning is the same as in "Step 3". In this embodiment... Pick The midpoint within. The angular deflection quantities of each order of the slowly varying kernel in this embodiment are as follows: Figure 5 , Figure 6 , Figure 7 , Figure 8 As shown.

[0063] The gradually varying core is not differentiable at some locations. This embodiment uses a numerical approximation solution obtained from discrete points in a sub-interval. It should be noted that the angular deflection variables of each order of the gradually varying core are determined by the seabed elastic parameters, and these variables are continuously updated and changed during the inversion process. Figure 5 , Figure 6 , Figure 7 , Figure 8 The seabed elastic parameters corresponding to the angular deflection quantities of the slowly varying core shown are as follows: (The text abruptly ends here, so the translation stops as well.) and high-speed seabed .

[0064] Step 5: Calculate the reflection coefficient vectors of each order of spherical wave. In this step... First-order reflection coefficient vector ,in and It is obtained from "step three" and "step four" respectively.

[0065] Step 6: Calculate the reflection coefficient of higher-order decoupled spherical waves. .

[0066] Figure 9 , Figure 10 The seabed elastic parameters used and Figure 5 , Figure 6 same, Figure 11 , Figure 12 The seabed elastic parameters used and Figure 7 , Figure 8 The results are identical. It is evident that the decoupled reflection coefficient and the accurate reflection coefficient of this invention are essentially in agreement, indicating that the proposed method for calculating the spherical wave reflection coefficient using oscillation-gradually varying core decoupling is effective. The calculation error of the decoupled reflection coefficient decreases with increasing decoupling order and with increasing number of sub-intervals. A balance between sub-interval discretization and calculation accuracy can be achieved by varying the decoupling order; calculation accuracy can be improved by increasing the decoupling order when the number of local sub-intervals is constant, or the number of sub-intervals can be reduced by increasing the decoupling order when calculation accuracy is constant. Especially for low-sound seabed, the effect of higher-order decoupling on suppressing calculation errors is particularly significant. Meanwhile, It is independent of the seabed elastic parameters, only This is related to seabed elastic parameters, which is the key to improving the efficiency of seabed elastic parameter inversion in the next step.

[0067] Step 7: Use the spherical wave reflection coefficients obtained from the higher-order decoupling in Step 6 to solve the subsequent inversion problem. For example... Figure 1 As shown in the shaded area, in the subsequent inversion problem solution process, it is only necessary to repeat "steps four, five, and six" after updating the seabed elastic parameters, without having to repeat the calculation of the local integral matrices of each order of oscillation kernel in step three.

[0068] This embodiment employs Markov chain Monte Carlo simulation to estimate the Bayesian posterior probability distribution of seabed elastic parameters. The inversion for each case utilizes two Markov chains computed in parallel. Each Markov chain iterates 12,000 times, retaining the stable data from the last 10,000 iterations, and merging the two Markov chains as the sampling result. The inversion under the oscillation-gradient kernel decoupling framework of this invention only requires updating the angular deflection parameters of each order of the gradient kernel during the iterative sampling calculation, eliminating the need to repeatedly calculate the local integral matrix of the oscillation kernel (see process below). Figure 1 (Shaded area). Traditional exact equation-driven inversion requires calculating the integral of the oscillation function in each iteration. Table 1 shows the computation time for the decoupled inversion of this invention and the traditional exact equation-driven inversion under different conditions.

[0069] Table 1. Comparison of computation time between the decoupled inversion of this invention and the traditional exact equation-driven inversion (unit: seconds)

[0070] ;

[0071] .

[0072] Table 1 shows the computation time of the above calculation process in the MATLAB R2023b environment. The time to the left of the slash " / " represents the time taken by an Intel(R) Core(TM) i9-10980XE CPU, and the time to the right of the slash " / " represents the time taken by a Hygon C86-3G (OPN:3350) CPU. The calculation of the local integral matrix of the oscillating kernel in the table was completed by 18 cores in parallel computation on the Intel CPU and by 8 cores in parallel computation on the Hygon CPU. Although the absolute time taken by different computing platforms will vary, the computation time of the method proposed in this invention is significantly less than that of the exact equation-driven inversion. The inversion results under different conditions are shown in Table 2. In the table, the first row of the cells for the seabed parameters is the expected posterior probability value, and the second row is the standard deviation of the absolute value. As can be seen from Table 2, the results obtained by the decoupled inversion method are basically consistent with the results of the traditional exact equation inversion.

[0073] Table 2 Comparison of Inversion Results

[0074] .

[0075] It is worth noting that although the inversion results in Table 2 are basically consistent under different conditions, when there is a significant error between the decoupled reflection coefficient and the accurate reflection coefficient (such as...), Figure 11 The posterior probability density distribution (as shown) will differ significantly. Errors in the decoupling reflection coefficient will lead to a difference in the density-P-wave velocity joint probability distribution (compared to...). Figures 13 to 15 (A) appears as a strip, with a joint probability distribution of shear wave velocity and longitudinal wave velocity (compare) Figures 13 to 15 (B) appears as a band, and the probability distribution of the longitudinal wave velocity shows a multi-peak state (compared to) Figures 13 to 15 (D); It can also be seen that the above phenomenon is weakened after the decoupling order changes from zero to second order. As the number of subintervals increases (from I=450 to I=1800 in this embodiment), the difference between the decoupled reflection coefficient and the exact reflection coefficient becomes smaller and smaller (D). Figure 12 When this difference decreases to a certain extent, its impact on the posterior probability distribution becomes negligible. Figure 16 , Figure 17 In this case, the slight difference in the posterior probability distribution only stems from the randomness of the sampling. Due to space limitations, the posterior probability distributions of the inversion results in other cases shown in Table 2 are not presented. In summary, the error in the decoupling reflection coefficients will cause changes in the shape of the posterior probability distribution, but it has little impact on its mean and variance. Increasing the number of subintervals or increasing the decoupling order can improve the accuracy of the decoupling reflection coefficients, thereby making the posterior probability distribution of the decoupling inversion results more consistent with the posterior probability distribution of the exact equation inversion results.

[0076] This invention decouples the oscillating kernel and the gradually varying kernel, avoiding redundant calculations of the oscillating kernel during the inversion process and significantly improving the calculation efficiency of the seabed spherical wave reflection coefficient during forward and inversion processes. Higher-order decoupling methods can achieve a balance between sub-interval discretization and computational accuracy by varying the decoupling order: computational accuracy can be improved by increasing the decoupling order when the number of local sub-intervals is constant, or the number of sub-intervals can be reduced when computational accuracy is constant.

[0077] Example 2

[0078] This embodiment also provides a high-order decoupling system for seismic spherical reflection wave oscillation-gradient kernel for seabed physical parameter inversion; the system includes a data input and initialization module, an integral interval decomposition module, a module for calculating the local integral matrix of each order oscillation kernel, a module for calculating the angular deflection vectors of each order of the gradually varying kernel, a module for calculating the spherical wave reflection coefficient vectors of each order, a module for calculating the spherical wave reflection coefficients of the high-order decoupled system, and a solution module for the subsequent inversion problem;

[0079] The data input module and initialization module execute the first step of the method; the integral interval decomposition module executes the second step; the module for calculating the local integral matrix of each oscillating kernel executes the third step; the module for calculating the angular deflection vectors of each order of the slowly varying kernel executes the fourth step; the module for calculating the spherical wave reflection coefficient vectors of each order executes the fifth step; the module for calculating the higher-order decoupled spherical wave reflection coefficient executes the sixth step; and the subsequent inversion problem solution module executes the sixth step. The specific steps are as follows:

[0080] The data input module and initialization module input observation data. Initialize seabed elastic parameters ;

[0081] The integration interval decomposition module decomposes the integration interval into several sub-intervals and sets the decoupling order;

[0082] The module for calculating the local integral matrix of each order of oscillation kernel calculates the local integral matrix of each order of oscillation kernel corresponding to the incident angle of the observed data.

[0083] The module for calculating the angular deflection parameters of each order of the slowly varying kernel calculates the angular deflection parameters of each order of the slowly varying kernel.

[0084] The module for calculating the reflection coefficient vector of each order of spherical wave calculates the reflection coefficient vector of each order of spherical wave; in this step, the first... First-order reflection coefficient vector ,in and The results are obtained by the module for calculating the local integral matrix of each order of oscillating kernel and the module for calculating the angular deflection vectors of each order of slowly varying kernel, respectively.

[0085] The module for calculating the higher-order decoupled spherical wave reflection coefficient calculates the higher-order decoupled spherical wave reflection coefficient. ;

[0086] The subsequent inversion problem solving module uses the higher-order decoupled spherical wave reflection coefficient in the higher-order decoupled spherical wave reflection coefficient calculation module to solve the subsequent inversion problem. In the process of solving the subsequent inversion problem, it is only necessary to repeatedly execute the modules for calculating the angle deflection vectors of the slowly varying kernel, the modules for calculating the spherical wave reflection coefficient vectors of each order, and the modules for calculating the spherical wave reflection coefficient vectors of each order after updating the seabed elastic parameters.

[0087] This embodiment also provides a high-order decoupling device for seismic spherical reflection wave oscillation-gradient kernel for seabed physical parameter inversion. The device is equipped with the above-mentioned decoupling system and executes the decoupling method described in Embodiment 1.

Claims

1. A method for decoupling spherical wave oscillations and slowly varying kernels for property inversion, characterized in that, The method steps are as follows: Step 1: Obtain observation data And input, initialize the seabed elastic parameters ; Step 2: Based on the observation data Determine the integration interval, decompose the integration interval into several sub-intervals, and set the decoupling order. Based on the observation data input in the first step Determine the integration interval Undetermined parameters in The integration interval is discretized into several sub-intervals, and then... To represent any one of the subintervals, ; Step 3: Calculate the local integral matrices of each order of oscillation kernel corresponding to the incident angle of the observed data; the local integral matrices of each order of oscillation kernel are used... It means that among them , representing zeroth order, first order, second order, ... Order; Local integral matrix of each order of oscillation kernel ;in, ... Indicates the angle of incidence corresponding to the observed data. ; Indicates a subinterval under a given incident angle. Internal oscillating nucleus of The value of a definite integral of order; ,in, The wave number of the seismic wave. Indicates the propagation distance of the wave field. It is a natural constant. express factorial, sub-interval Taylor series expansion points within; Step 4: Calculate the angular deflection coefficients of each order of the slowly varying kernel; sequentially for all subintervals in Step 2. Internally determined by the current seabed elastic parameters Calculate the plane wave reflection coefficient For the angle of incidence The partial derivatives of each order are used to obtain the angular partial derivatives of the slowly varying kernel. ; Step 5: Calculate the reflection coefficient vectors of each order of spherical wave; in this step, the first... First-order reflection coefficient vector ,in, and Obtained from steps three and four respectively; Step 6: Calculate the reflection coefficient of higher-order decoupled spherical waves. ; Step 7: Use the higher-order decoupled spherical wave reflection coefficients obtained in Step 6 to solve the subsequent inversion problem.

2. The spherical wave oscillation-gradiently varying kernel high-order decoupling method for property inversion according to claim 1, characterized in that, The value of does not need to be strictly calculated; it is only necessary to ensure that the integral outside the interval is much smaller than the integral inside the interval; the union of all subintervals forms the integration interval. The intersection of any two subintervals is an empty set. Let the decoupling order be denoted as . .

3. The spherical wave oscillation-gradiently varying kernel high-order decoupling method for property inversion according to claim 1, characterized in that, In the fourth step, when calculating the angular deflection quantities of each order of the slowly varying core for the first time, the current seabed elastic parameters are... Refers to the seabed elastic parameters initialized in the first step. In subsequent calculations of the angular deflection derivatives of the slowly varying core, the current seabed elastic parameters are determined during the iterative solution process. Refers to the seabed elasticity parameters after the seventh step update. ; , , , representing zeroth order, first order, second order, ... Rank.

4. The spherical wave oscillation-gradiently varying kernel high-order decoupling method for property inversion according to claim 1, characterized in that, In the seventh step, in the subsequent inversion problem solution process, it is only necessary to repeat the fourth, fifth and sixth steps after updating the seabed elastic parameters, without having to repeat the calculation of the local integral matrix of each order of oscillation kernel in the third step, thus avoiding repeated integration calculation of the oscillation function.

5. A high-order decoupling device for spherical wave oscillation-gradiently varying kernel for property inversion, characterized in that, The device is equipped with a spherical wave oscillation-slowly varying kernel high-order decoupling system for property inversion, and the decoupling system performs the decoupling method according to any one of claims 1-4; The decoupling system includes a data input and initialization module, an integral interval decomposition module, a module for calculating the local integral matrix of each order of oscillating kernel, a module for calculating the angular deflection vectors of each order of slowly varying kernel, a module for calculating the spherical wave reflection coefficient vectors of each order, a module for calculating the spherical wave reflection coefficients of higher-order decoupled systems, and a solution module for the subsequent inversion problem. The data input and initialization module is used to input observation data. Initialize seabed elastic parameters ; The integral interval decomposition module is used to decompose the integral interval into several sub-intervals and set the decoupling order; The module for calculating the local integral matrix of each oscillation kernel is used to calculate the local integral matrix of each oscillation kernel corresponding to the incident angle of the observed data. The module for calculating the angular deflection parameters of each order of the slowly varying kernel is used to calculate the angular deflection parameters of each order of the slowly varying kernel. The module for calculating the reflection coefficient vector of spherical waves of each order is used to calculate the reflection coefficient vector of spherical waves of each order. First-order reflection coefficient vector ,in and The results are obtained by the module for calculating the local integral matrix of each order of oscillating kernel and the module for calculating the angular deflection vectors of each order of slowly varying kernel, respectively. The module for calculating the higher-order decoupled spherical wave reflection coefficient is used to calculate the higher-order decoupled spherical wave reflection coefficient. ; The module for solving the subsequent inversion problem is used to solve the subsequent inversion problem by utilizing the higher-order decoupled spherical wave reflection coefficients from the module for calculating higher-order decoupled spherical wave reflection coefficients.

Images