Complex thin interbedded inversion method for single thin layer of VTI and related equipment

Abstract

The application provides a complex thin interbed inversion method for a VTI single thin layer, comprising the following steps: constructing an initial model parameter matrix according to the parameterization of the underground medium; calculating a reflection coefficient matrix of the underground medium according to the initial model parameter matrix, and using the reflection coefficient matrix to forward simulate a synthetic seismic record. The synthetic seismic record comprises a PP wave synthetic record and a PS wave synthetic record. An inversion target function is constructed according to the actual seismic record and the synthetic seismic record of the target underground medium, and the model parameter matrix and the synthetic seismic record are used in the construction of the target function. The actual seismic record comprises a PP wave actual record and a PS wave actual record. The inversion target function is solved, the initial parameter model is updated and iterated, and thin interbed parameter information is obtained. Thus, the joint inversion of the PP wave and the PS wave of the complex thin interbed with the VTI single thin layer is realized, and the inversion precision of the complex thin interbed with the VTI single thin layer is improved.

Description

Technical Field

[0001] This application relates to the field of seismic exploration technology, and in particular to a thin interlayer inversion method. Background Technology

[0002] In recent years, with the continuous deepening of oil and gas exploration, the wave field propagation laws and reservoir property parameter inversion in complex structures have been studied more extensively and in-depth. Against this research background, the refined exploration and efficient development of thin interbedded reservoirs has become one of the important goals of oil and gas exploration and development, and the accuracy requirements for characterizing their property parameters have also increased accordingly. Traditional seismic inversion techniques face severe challenges when dealing with thin interbedded reservoirs: on the one hand, the thickness of a single layer in thin interbedded reservoirs is generally lower than the conventional seismic resolution limit (usually less than 1 / 4 wavelength); on the other hand, the complex layered sedimentary structure causes transmission loss, destructive interference, and dispersion attenuation of seismic wave energy, leading to insufficient vertical resolution and reduced accuracy of property parameter inversion results.

[0003] Traditional AVO inversion relies on the Zoeppritz equation and its approximations to calculate various parameters of the subsurface medium, and this method is widely used. However, the Zoeppritz equation is based on the single-interface assumption, and when dealing with the inversion of thin interlayered structures, it cannot directly calculate the reflection coefficient of each interface within the thin interlayered structure. This makes it difficult to accurately characterize the complex internal structure and wave field propagation laws of the thin interlayered structure, which may lead to calculation errors, especially in regions with complex thin interlayered structures.

[0004] Other studies utilize recursive matrix methods for inversion of complex media and thin interbedded layers. Compared to inversion methods based on the Zoeppritz equation, these methods can more accurately calculate complex wavefield characteristics such as transmission loss, interlayer multiples, and wave mode conversions, making them more suitable for inverting thin interbedded layers. There are many studies on AVO methods based on recursive matrices, with the most widely used being the reflectance forward modeling (KRM) method proposed by Kennett (1983) using the Kennett recursive matrix. Most existing studies treat thin interbedded layers as multilayered isotropic media. However, in some complex thin interbedded layers, individual thin layers exhibit short-term cyclic characteristics, and the grain size within these layers shows rhythmicity and periodicity. Some thin layers composed of multiple extremely thin layers also exhibit regularly varying lithologies. These individual thin layers themselves possess VTI characteristics; treating them as isotropic layers fails to accurately characterize the physical parameters of the thin interbedded layers. Summary of the Invention

[0005] In view of this, the purpose of this application is to provide a thin interlayer inversion method, apparatus, electronic device and storage medium.

[0006] To achieve the above objectives, this application provides a thin interlayer inversion method, comprising:

[0007] An initial model parameter matrix is ​​constructed based on the parameterization of the subsurface medium. The reflection coefficient matrix of the subsurface medium is calculated based on this initial model parameter matrix, and the synthetic seismic record is simulated using this reflection coefficient matrix in forward modeling. The synthetic seismic record includes PP wave and PS wave synthetic records. Based on the actual seismic record of the target subsurface medium and the synthetic record, the corresponding forward modeling operator and algorithm (a second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm) are selected. The actual seismic record includes actual PP wave and PS wave records. The inversion objective function is solved, and the model parameters are iteratively updated. The updated parameter model is used as the new initial model parameters, and the above steps are repeated until the thin interbedded layer parameter information is obtained. This achieves the joint inversion of PP and PS waves in complex thin interbedded layers with a single VTI layer, improving the inversion accuracy for such layers.

[0008] In some embodiments, the reflection coefficient matrix of the subsurface medium is calculated based on the initial model parameters, and the reflection coefficient matrix is ​​used to perform forward modeling to synthesize seismic records. Specifically, this includes: parameterizing the subsurface medium to obtain the initial model parameter matrix: m = [V P0 V S0 ρδε] T , where V P0 V S0 ρ, δ, and ε represent the P-wave velocity and S-wave velocity density parameter vectors and two Thomson parameter vectors, respectively. The reflection coefficient matrix is ​​calculated based on the model parameters. The synthetic seismic record is then calculated using the reflection coefficient matrix and a preset wavelet matrix. The synthetic seismic record is calculated using the following formula: d(m,θ,t)=W(t)R(m,θ,t); where d(m,θ,t) represents the synthetic seismic record, W(t) represents the wavelet matrix, R(m,θ,t) represents the reflection coefficient matrix, m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0009] In some embodiments, the calculation of the reflection coefficient matrix specifically includes: calculating the PP wave reflection coefficient and PS wave reflection coefficient corresponding to the model parameters using the Kennett second-order approximation formula based on the initial model parameter matrix, thereby obtaining the reflection coefficient matrix; the slow-frequency domain second-order approximation reflection coefficient is calculated by the following formula: Where ω represents angular frequency and p represents horizontal slowness. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n-1 for a downlink wave in the VTI medium. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for a downlink wave in the VTI medium. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the downlink wave. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for the upward wave. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the upward wave. Let be the phase shift factor as the wave propagates from interface n-1 to interface n. When calculating synthetic seismic records, a second-order approximate reflection coefficient in the angle-time domain is required, which is calculated using the following formula: m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0010] In some embodiments, constructing the inversion objective function based on the synthetic seismic record specifically includes: calculating the difference between the synthetic seismic record and the measured data to obtain a biased seismic record; calculating the difference between the model parameter matrix and the measured data to obtain biased model parameters; and constructing the inversion objective function based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm, wherein the model parameter matrix and the synthetic seismic record are used in the construction of the objective function; the inversion objective function is expressed by the following formula: J(m)=||(H+λ1I)Δm-J T d ref || 2 +λ2||Δm+m ref ||,; where λ1 represents the damping coefficient, λ2 represents the sparse constraint coefficient, J represents the Jacobian matrix, and H represents the Hessian matrix H=J T J and I denote the identity matrix, d ref The deviation seismic record represents the value of the seismic data, and Δm represents the update amount of the model parameter matrix. ref This represents the parameters of the deviation model. The J-expansion form under multi-wave joint inversion is:

[0011] In some embodiments, solving the inversion objective function and updating the model parameter matrix to obtain the complex thin interlayer parameter information sandwiched with a single VTI thin layer specifically includes: setting the update amount of the model parameters; iterating the inversion objective function according to the update amount of the model parameters; and, in response to the inversion objective function approaching 0, using the current model parameters as the thin interlayer parameter information; the update amount of the model parameters is expressed by the following formula: Where, d PP This represents the PP wave synthesis record, d PS This represents the PS wave synthesis record, s PPThis represents the measured data of the PP wave, s PS This represents the measured data of the PS wave.

[0012] A thin interlayer inversion device includes: a parameter construction module for constructing an initial model parameter matrix based on preset seismic parameters; a seismic calculation module for calculating a reflection coefficient matrix based on the initial model parameter matrix, and then calculating a synthetic seismic record; wherein the synthetic seismic record includes a PP wave synthetic record and a PS wave synthetic record; an inversion module for constructing an inversion objective function based on a second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm, wherein the objective function is constructed using the model parameter matrix and the synthetic seismic record; and a solution module for solving the inversion objective function and updating the model parameter matrix to obtain the parameter information of the complex thin interlayer sandwiched with a single VTI thin layer.

[0013] In some embodiments, the seismic calculation module specifically includes: a forward modeling module, which calculates the PP wave reflection coefficient and PS wave reflection coefficient matrix corresponding to the model parameters according to Kennett's second-order approximation formula; and calculates the PP wave composite record and the PS wave composite record using the reflection coefficient matrix and a preset wavelet matrix; the composite seismic record is calculated by the following formula: d(m,θ,t)=W(t)R(m,θ,t); where d(m,θ,t) represents the composite seismic record, W(t) represents the wavelet matrix, R(m,θ,t) represents the reflection coefficient matrix, m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0014] An electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the method as described in any of the preceding descriptions.

[0015] A non-transitory computer-readable storage medium storing computer instructions for causing a computer to perform any of the methods described above.

[0016] As can be seen from the above, the inversion method for complex thin interbedded layers with VTI single thin layers provided in this application constructs an initial model parameter matrix, calculates the corresponding PP wave reflection coefficient and PS wave reflection coefficient matrices based on the initial model parameter matrix, calculates the PP wave composite record and the PS wave composite record using the reflection coefficient matrix and a preset wavelet matrix, and constructs an inversion objective function based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm. The objective function is constructed using the model parameter matrix and the composite seismic record. The objective function is solved and the model parameters are iterated, thereby achieving joint inversion of PP waves and PS waves and improving the inversion accuracy for complex thin interbedded layers with VTI single thin layers. Attached Figure Description

[0017] To more clearly illustrate the technical solutions in this application or related technologies, the drawings used in the description of the embodiments or related technologies will be briefly introduced below. Obviously, the drawings described below are only embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0018] Figure 1 A schematic flowchart of the thin interlayer inversion method provided in this application embodiment;

[0019] Figure 2 This is a schematic diagram of an AVA gather used to simulate calculation results in an embodiment of this application;

[0020] Figure 3 The inversion result obtained according to the method of this application is shown in the figure.

[0021] Figure 4 The image shows the inversion results obtained from the exact Zoeppritz equation.

[0022] Figure 5 This is a schematic diagram of a noisy AVA gather used to simulate calculation results in an embodiment of this application;

[0023] Figure 6 This is an inversion result diagram obtained by inverting an AVA gather containing noise according to the method of this application;

[0024] Figure 7 This is a schematic diagram of the thin interlayer inversion device provided in the embodiments of this application;

[0025] Figure 8 This is a schematic diagram of a more specific electronic device hardware structure provided in this application. Detailed Implementation

[0026] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with specific embodiments and the accompanying drawings.

[0027] It should be noted that, unless otherwise defined, the technical or scientific terms used in the embodiments of this application should have the ordinary meaning understood by one of ordinary skill in the art to which this application pertains. The terms "first," "second," and similar terms used in the embodiments of this application do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed after the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are only used to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0028] It is understood that before using the technical solutions of the various embodiments in this disclosure, users will be informed of the type, scope of use, and usage scenarios of the personal information involved in an appropriate manner, and user authorization will be obtained.

[0029] For example, upon receiving a user's active request, a prompt message is sent to the user to explicitly inform them that the requested operation will require the acquisition and use of the user's personal information. This allows the user to independently choose, based on the prompt message, whether to provide personal information to the software or hardware such as electronic devices, applications, servers, or storage media performing the operations of this disclosed technical solution.

[0030] As an optional but not limited implementation, in response to a user's active request, sending a prompt message to the user can be done via a pop-up window, where the prompt message can be presented in text format. Furthermore, the pop-up window can also include a selection control allowing the user to choose "agree" or "disagree" to provide personal information to the electronic device.

[0031] It is understood that the above notification and user authorization process are merely illustrative and do not constitute a limitation on the implementation of this disclosure. Other methods that comply with relevant laws and regulations may also be applied to the implementation of this disclosure.

[0032] Thomsen parameters are a set of parameters in seismic geophysics used to describe the elastic properties of rocks. These parameters describe the relationship between elastic wave velocity and density in isotropic rocks. Thomsen parameters include ε (epsilon), δ (delta), and γ (gamma), which describe the elastic anisotropy, inelastic effects, and isotropic properties of rocks, respectively. This application uses only the ε and δ parameters.

[0033] PP wave: This refers to the P-wave (compression wave) emitted from the seismic source that reflects back to the Earth's surface after encountering a subsurface interface. When a P-wave encounters a subsurface interface, some of its energy is reflected back to the surface, forming a PP wave. PP waves are mainly used to detect the density and elastic properties of subsurface rock layers and are one of the most common waveforms in seismic exploration.

[0034] PS waves are waves that, after being reflected from subsurface interfaces, are converted into S-waves (shear waves) and return to the surface. S-waves propagate in solids but not in fluids, thus being more sensitive to the internal structure of rock fabrics. PS waves, after reflection from subsurface interfaces, provide crucial information about the physical properties of rocks. Compared to PP waves, PS waves offer higher resolution and accuracy in imaging fluid-saturated rocks.

[0035] Ricker wavelet matrix: It is a signal with a definite start time, finite energy and a certain duration, and it is the basic unit in seismic records.

[0036] like Figure 1 As shown, this application provides a thin interlayer inversion method, including:

[0037] Step S1: Construct an initial model parameter matrix based on preset seismic parameters.

[0038] In this embodiment, the model parameter matrix includes P-wave velocity, S-wave velocity, density parameter vector, and Thomson parameter vector. In a thin, anisotropic layer, the model parameter matrix is ​​represented in the following format:

[0039] m = [V P0 V S0 ρδε]T;

[0040] Among them, V P0 V represents the longitudinal wave velocity. S0 ρ represents the transverse wave velocity, δ represents the density parameter vector, δ represents the elastic anisotropy in the Thomsen parameters, and ε represents the inelastic effects in the Thomsen parameters.

[0041] Step S2: Calculate the reflection coefficient matrix based on the initial model parameter matrix; calculate the synthesized seismic record based on the wave reflection coefficient matrix and the preset wavelet matrix. The synthesized seismic record includes PP wave synthesized records and PS wave synthesized records.

[0042] Step S3: Based on the actual seismic record of the target subsurface medium and the synthetic record, select the corresponding forward modeling operator and algorithm (the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm) to construct the inversion objective function. The objective function is constructed using the model parameter matrix and the synthetic seismic record.

[0043] Step S4: Solve the inversion objective function to obtain the thin interlayer parameter information.

[0044] The inversion objective function is the optimization objective for inverting seismic waves. The parameters in the inversion process are adjusted according to the inversion objective function so that the thin interbedded parameter information obtained by inversion gradually approaches the actual geological conditions.

[0045] In some embodiments, step S2 specifically includes:

[0046] Step S21: Calculate the PP wave reflection coefficient and PS wave reflection coefficient matrix corresponding to the model parameters using the Kennett second-order approximation formula based on the initial model parameter matrix.

[0047] Step S22: Calculate the synthetic seismic record based on the reflection coefficient matrix and the preset wavelet matrix.

[0048] In this embodiment, the wavelet matrix is ​​the Ricker wavelet matrix.

[0049] Synthetic seismic records are calculated using the following formula:

[0050] d(m,θ,t)=W(t)R(m,θ,t).

[0051] Where d(m,θ,t) represents the synthetic seismic record, W(t) represents the wavelet matrix, R(m,θ,t) represents the reflection coefficient matrix, m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0052] In some embodiments, step S21 specifically includes:

[0053] Step S211: Calculate the PP wave reflection coefficient and PS wave reflection coefficient corresponding to the model parameters according to Kennett's second-order approximation formula.

[0054] In this embodiment, the second-order approximation of the Kennett reflection coefficient equation, proposed by Huang and Lu (2022), is used in the inversion process for calculating the reflection coefficient of a VTI thin layer. By transforming the isotropic single-interface reflection coefficient term in the second-order approximation of the Kennett equation into a single-interface reflection coefficient term under the VTI medium assumption, a second-order approximation of the Kennett equation applicable to VTI media is obtained. This method is suitable for calculating the reflection coefficient matrix of a VTI single thin layer with PP and PS waves. Then, using the recursive method proposed by Kennett (1983), the second-order approximate reflection coefficient of any interface n of the thin interlayer system with a VTI single thin layer sandwiched in the slow-frequency domain is derived, where the incident wave is a P-wave or an S-wave.

[0055] Step S212: Calculate the synthesized seismic record based on the wave reflection coefficient and the preset wavelet matrix.

[0056] In this embodiment, the reflection coefficient matrix is ​​represented by the following formula:

[0057]

[0058] Where m represents the model parameter matrix, ω represents the angular frequency, and p represents the horizontal slowness. Indicates the reflection coefficient of PP waves. Represents the PS wave reflection coefficient. Indicates the reflection coefficient of SP waves. This represents the reflection coefficient of SS waves.

[0059] Preferably, in the inversion process, since only the reflection coefficients of the PP wave and PS wave need to be calculated, the above equation can be simplified in practical applications to obtain:

[0060]

[0061] In some embodiments, step S211 specifically includes:

[0062] Step S201: Construct the second-order approximate reflection coefficient matrix according to Kennett's second-order approximation formula. The second-order approximate reflection coefficients are calculated by the following formula:

[0063]

[0064] Where ω represents angular frequency and p represents horizontal slowness. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n-1 for a downlink wave in the VTI medium. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for a downlink wave in the VTI medium. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the downlink wave. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for the upward wave. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the upward wave. Let be the phase shift factor as the wave propagates from interface n-1 to interface n.

[0065] in, Used to represent the first-order reflection coefficient of interface n. It includes the primary reflection coefficient of interface n-1 and the first-order and second-order multiple reflection coefficients of layer n-1.

[0066] Under the isotropic assumption, the single-interface coefficients in Kennett's second-order approximation are displacement anti-transmission coefficients and transmission coefficients. Huang and Lu (2020), under the VTI medium assumption, replaced the isotropic single-interface anti-transmission coefficients and transmission coefficients with those calculated based on Graebner's formula (1992). It is important to note that the single-interface coefficients in Graebner's formula are displacement anti-transmission coefficients and transmission coefficients, which differ from those in Kennett's formula. To achieve the correct substitution between the two types of coefficients, the anti-transmission coefficients calculated by Graebner must be multiplied by the ratio of the corresponding wave velocities on both sides of the interface, thus converting them into a displacement anti-transmission coefficient matrix that conforms to the VTI medium assumption. and

[0067]

[0068] Where the element r in the matrix ** t ** V represents the single-interface reflectance and transmittance coefficients for different modes calculated using the Graebner formula. P V S Represents the longitudinal and transverse wave velocities. The superscripts n and n-1 of the element represent the nth and n-1th media layers, respectively.

[0069] In this embodiment, the phase shift factor is used to ensure the phase continuity of the wave as it passes through the interface in the propagation matrix method. The phase shift factor matrix in the VTI medium is expanded as follows:

[0070]

[0071] q P (θ P ) and q S (θ S ) represents the vertical slowness of the P-wave and S-wave, and h(n-1) represents the thickness of layer n-1. For q P (θ P ) and qS (θ S For all of them, there are:

[0072]

[0073] in:

[0074]

[0075] c 11 c 33 c 55 c 13 It is an independent stiffness coefficient (Thomsen, 1986):

[0076]

[0077] c 11 =2c 33 ε+c 33 ;

[0078]

[0079] For weakly anisotropic media, according to Christoffel's equations, the ray parameters of the P-wave are:

[0080]

[0081] The ray parameters of the S-wave are:

[0082]

[0083] Since the reflection coefficient matrix calculated using the second-order approximate reflection coefficient formula is in the slowness-frequency domain, while the reflection coefficients used for inversion in this application should be in the angle-time domain, it is necessary to transform the second-order approximate reflection coefficient calculation formula. According to the inverse Fourier transform and Snell's law, the formula for the total reflection coefficient in the angle-time domain can be obtained as follows:

[0084]

[0085] Step S203: Represent the preset PS wave AVA gather in the PP wave time domain to obtain the converted PS wave AVA gather.

[0086] Preferably, the PS AVA gather can be compressed to the PP time domain using a conversion factor, and the corresponding conversion factor for each sampling point is:

[0087] Where i represents different incident angles, and j represents different sampling times. V P,ij V represents the P-wave velocity sampled at incident angle i and time j;S,ij This represents the S-wave velocity at incident angle i and time j.

[0088] In some embodiments, step S3 specifically includes:

[0089] Step S31: Calculate the difference between the synthetic seismic record and the measured data to obtain the biased seismic record.

[0090] In this embodiment, the deviation seismic record is calculated according to the following formula:

[0091]

[0092] Where, d PP This indicates the PP wave synthesis record, d PS Indicates PS wave synthesis record, s PP This represents the actual seismic data of PP waves, s PS This represents actual seismic data from PS waves.

[0093] Step S32: Calculate the difference between the model parameter matrix and the measured data to obtain the deviation model parameters.

[0094] The deviation model parameters are calculated according to the following formula:

[0095] m ref =m0-m true ;

[0096] Where m0 represents the model parameter matrix currently involved in the calculation; m true This represents the model parameter matrix obtained from actual measurements, typically using well logging data.

[0097] Step S33: Construct the inversion objective function based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm. The objective function is constructed using the model parameter matrix and the synthetic seismic record.

[0098] The inversion objective function is expressed by the following formula:

[0099] J(m)=||(H+λ1I)Δm-J T d ref || 2 +λ2||Δm+m ref ||.

[0100] Where λ1 represents the damping coefficient, λ2 represents the sparse constraint coefficient, J represents the Jacobian matrix, and H represents the Hessian matrix H=J T J and I denote the identity matrix, d ref The variable represents the biased seismic record, and Δm represents the update amount of the model parameter matrix.ref This represents the parameters of the deviation model.

[0101] The inversion objective function is constructed based on the LM algorithm under the ideal condition that the forward seismic record has zero noise.

[0102] The preferred form of J-expansion under multi-wave joint inversion is:

[0103]

[0104] In some embodiments, step S4 specifically includes:

[0105] Step S41: Set the model parameter update amount and iterate the inversion objective function based on the model parameter update amount.

[0106] Step S42: In response to the inversion objective function approaching 0, the current model parameters are used as thin interlayer parameter information.

[0107] In this embodiment, when the inversion objective function approaches 0, the model parameter matrix m0 currently involved in the calculation and the measured model parameter matrix m0 are compared. true Extremely close, that is, at this moment:

[0108]

[0109] As an optional implementation, an inversion threshold can be set for the inversion objective function. When the value of the inversion objective function is lower than the inversion threshold, the inversion is stopped, thereby improving the inversion efficiency.

[0110] The model parameter update amount is expressed by the following formula:

[0111]

[0112] Where, d PP This indicates the PP wave synthesis record, d PS Indicates PS wave synthesis record, s PP This represents the measured data of the PP wave, s PS This represents the measured data of the PS wave.

[0113] As can be seen from the above embodiments of this application, the thin interlayer inversion method provided by this application constructs an initial model parameter matrix, then calculates synthetic seismic records including PP wave synthetic records and PS wave synthetic records based on the initial model parameter matrix, and constructs an inversion objective function based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm. The objective function is constructed using the model parameter matrix and synthetic seismic records. The objective function is solved and the model parameter matrix is ​​updated, thereby realizing the joint inversion of PP waves and PS waves and improving the inversion accuracy of complex thin interlayers with VTI single thin layers.

[0114] As an exemplary embodiment, this application uses a thin interlayer model with a single VTI layer sandwiched between it to verify the inversion method. The parameters are shown in the table below:

[0115]

[0116] Forward modeling was performed on the model, with a frequency calculation range of 0–105 Hz and a source wavelet of 30 Hz Ricker wavelet. The AVA gathers of the PP and PS waves of the model were obtained, as follows: Figure 2 As shown, considering that the incident angle range of actual data angle gathers is usually narrow, the angle range of the angle gathers in this thin interlayer model is 5° to 30°, with an interval of 5°.

[0117] Smooth the real model to obtain the inversion initial model. Figure 3 The multi-wave AVA gathers shown are the target data for inversion. During inversion, the reflection coefficients were calculated using both the Kennett second-order approximation formula applicable to VTI media and the exact Zoeppritz equation for comparison.

[0118] in, Figure 3 This indicates the inversion result obtained by performing the inversion according to the method described in this application. Figure 3 (a) represents the longitudinal wave velocity. Figure 3 (b) represents the transverse wave velocity. Figure 3 (c) represents density. Figure 3 (d) represents the value of δ. Figure 3 (e) represents the ε value. The black line represents the initial model, the red line represents the true model, and the blue line represents the inversion result. It can be seen that the inversion result obtained by the method of this application has a high degree of fit with the true value.

[0119] in, Figure 4 The results of the inversion based on the exact Zoeppritz equation are as follows: Figure 4 (a) represents the longitudinal wave velocity. Figure 4 (b) represents the transverse wave velocity. Figure 4(c) represents density. The black line represents the initial model, the red line represents the true model, and the blue line represents the inversion result. It can be seen that the inversion result based on the Zoeppritz equation will have some fluctuations. This may be due to the complex structure inside the thin interlayer and the influence of VTI characteristics on the wave field. The Zoeppritz equation does not take their influence into account when calculating the reflection coefficient, so the accuracy is reduced.

[0120] As an exemplary embodiment, to test the noise resistance of the inversion method, this application also designs a thin interlayer model synthesis record with added random noise signal, such as... Figure 5 As shown.

[0121] Using a smoothing model as the initial model, inversion is performed based on the method provided in this application.

[0122] Figure 6 As shown in the inversion results, the results fluctuate due to random noise. However, at the target layer, the inversion results still reflect the true parameters of the formation well, indicating that the inversion method has good noise resistance.

[0123] in, Figure 6 This represents the inversion result obtained by using the method provided in this application to invert noisy data: Figure 6 (a) represents the longitudinal wave velocity. Figure 6 (b) represents the transverse wave velocity. Figure 6 (c) represents density. Figure 6 (d) represents the value of δ. Figure 6 (e) represents the ε value. The black line represents the initial model, the red line represents the true model, and the blue line represents the inversion result. This result demonstrates that the method proposed in this application is more capable of inverting the complex internal wavefields of thin, interlayered structures, and can produce more accurate physical property parameter models.

[0124] It should be noted that the method in this embodiment can be executed by a single device, such as a computer or server. The method can also be applied in a distributed scenario, where multiple devices cooperate to complete the task. In such a distributed scenario, one of these devices may execute only one or more steps of the method in this embodiment, and the multiple devices will interact with each other to complete the method described.

[0125] It should be noted that the above description describes some embodiments of this application. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in a different order than that shown in the above embodiments and still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0126] Based on the same inventive concept, corresponding to any of the above embodiments, this application also provides a thin interlayer inversion device.

[0127] refer to Figure 7 Thin interlayer inversion device, comprising:

[0128] The parameter construction module 100 is used to construct the model parameter matrix based on preset seismic parameters.

[0129] The seismic calculation module 200 calculates the reflection coefficient matrix of the subsurface medium based on the initial model parameter matrix, and uses the reflection coefficient matrix to perform forward modeling to synthesize seismic records; wherein, the synthesized seismic records include PP wave synthesized records and PS wave synthesized records.

[0130] Inversion module 300 constructs an inversion objective function based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm. The objective function is constructed using the model parameter matrix and the synthetic seismic record.

[0131] The solver module 400 is used to solve the inversion objective function to obtain the thin interlayer parameter information.

[0132] In some embodiments, the seismic calculation module 200 specifically includes:

[0133] Forward modeling module 210 calculates the PP wave reflection coefficient and PS wave reflection coefficient matrix corresponding to the model parameters according to Kennett's second-order approximation formula;

[0134] The PP wave synthesis record and the PS wave synthesis record are obtained by calculation using the reflection coefficient matrix and the preset wavelet matrix.

[0135] Synthetic seismic records are calculated using the following formula:

[0136] d(m,θ,t)=W(t)R(m,θ,t).

[0137] Where d(m,θ,t) represents the synthetic seismic record, W(t) represents the wavelet matrix, R(m,θ,t) represents the reflection coefficient matrix, m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0138] For ease of description, the above devices are described in terms of function, divided into various modules. Of course, in implementing this application, the functions of each module can be implemented in one or more software and / or hardware.

[0139] The apparatus of the above embodiments is used to implement the corresponding thin interlayer inversion method in any of the foregoing embodiments, and has the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0140] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this application also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the thin cross-layer inversion method described in any of the above embodiments.

[0141] Figure 8 This embodiment illustrates a more specific hardware structure of an electronic device, which may include a processor 1010, a memory 1020, an input / output interface 1030, a communication interface 1040, and a bus 1050. The processor 1010, memory 1020, input / output interface 1030, and communication interface 1040 are interconnected internally via the bus 1050.

[0142] The processor 1010 can be implemented using a general-purpose CPU (Central Processing Unit), microprocessor, application-specific integrated circuit (ASIC), or one or more integrated circuits, and is used to execute relevant programs to implement the technical solutions provided in the embodiments of this specification.

[0143] The memory 1020 can be implemented in the form of ROM (Read Only Memory), RAM (Random Access Memory), static storage device, dynamic storage device, etc. The memory 1020 can store the operating system and other applications. When the technical solutions provided in the embodiments of this specification are implemented by software or firmware, the relevant program code is stored in the memory 1020 and is called and executed by the processor 1010.

[0144] The input / output interface 1030 is used to connect input / output modules to realize information input and output. Input / output modules can be configured as components within the device (not shown in the figure) or externally connected to the device to provide corresponding functions. Input devices may include keyboards, mice, touchscreens, microphones, various sensors, etc., while output devices may include displays, speakers, vibrators, indicator lights, etc.

[0145] The communication interface 1040 is used to connect a communication module (not shown in the figure) to enable communication between this device and other devices. The communication module can communicate via wired means (such as USB, Ethernet cable, etc.) or wireless means (such as mobile network, WIFI, Bluetooth, etc.).

[0146] Bus 1050 includes a pathway for transmitting information between various components of the device, such as processor 1010, memory 1020, input / output interface 1030, and communication interface 1040.

[0147] It should be noted that although the above-described device only shows the processor 1010, memory 1020, input / output interface 1030, communication interface 1040, and bus 1050, in specific implementations, the device may also include other components necessary for normal operation. Furthermore, those skilled in the art will understand that the above-described device may only include the components necessary for implementing the embodiments of this specification, and not necessarily all the components shown in the figures.

[0148] The electronic devices described above are used to implement the corresponding thin interlayer inversion method in any of the foregoing embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0149] Based on the same inventive concept, corresponding to the methods of any of the above embodiments, this application also provides a non-transitory computer-readable storage medium storing computer instructions for causing the computer to execute the thin cross-layer inversion method as described in any of the above embodiments.

[0150] The computer-readable medium of this embodiment includes permanent and non-permanent, removable and non-removable media, and information storage can be implemented by any method or technology. Information can be computer-readable instructions, data structures, program modules, or other data. Examples of computer storage media include, but are not limited to, phase-change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, CD-ROM, digital versatile optical disc (DVD) or other optical storage, magnetic tape, magnetic magnetic disk storage or other magnetic storage devices, or any other non-transfer medium that can be used to store information accessible by a computing device.

[0151] The computer instructions stored in the storage medium of the above embodiments are used to cause the computer to execute the thin interlayer inversion method as described in any of the above embodiments, and have the beneficial effects of the corresponding method embodiments, which will not be repeated here.

[0152] It should be noted that the embodiments of this application can also be further described in the following ways:

[0153] A thin interlayer inversion method includes:

[0154] The model parameter matrix is ​​constructed based on the preset earthquake parameters.

[0155] The reflection coefficient matrix is ​​calculated based on the model parameter matrix, and then the synthetic seismic record is calculated based on the reflection coefficient matrix and the preset wavelet matrix. The synthetic seismic record includes PP wave synthetic records and PS wave synthetic records.

[0156] The inversion objective function is constructed based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm, utilizing the model parameter matrix and synthetic seismic records.

[0157] The inversion objective function is solved and the model parameter matrix is ​​updated to obtain the thin interlayer parameter information.

[0158] In some embodiments, the calculation of synthetic seismic records based on the model parameter matrix specifically includes:

[0159] The reflection coefficient matrix is ​​calculated based on the model parameters.

[0160] The synthetic seismic record is calculated based on the reflection coefficient matrix and the preset wavelet matrix.

[0161] Synthetic seismic records are calculated using the following formula:

[0162] d(m,θ,t)=W(t)R(m,θ,t).

[0163] Where d(m,θ,t) represents the synthetic seismic record, W(t) represents the wavelet matrix, R(m,θ,t) represents the reflection coefficient matrix, m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0164] In some embodiments, forward modeling is performed on the PP wave composite record and the PS wave composite record, specifically including:

[0165] The PP wave reflection coefficient and PS wave reflection coefficient corresponding to the model parameters are calculated based on Kennett's second-order approximation formula.

[0166] The seismic record is calculated and synthesized based on the wave reflection coefficient and the preset wavelet matrix.

[0167] In some embodiments, the calculation of the PP wave reflection coefficient and PS wave reflection coefficient corresponding to the model parameters according to the Kennett second-order approximation formula specifically includes:

[0168] Construct a second-order approximate reflection coefficient matrix based on Kennett's second-order approximation formula.

[0169] The second-order approximate reflection coefficient is calculated using the following formula:

[0170]

[0171] Where ω represents angular frequency and p represents horizontal slowness. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n-1 for a downlink wave in the VTI medium. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for a downlink wave in the VTI medium. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the downlink wave. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for the upward wave. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the upward wave. Let be the phase shift factor as the wave propagates from interface n-1 to interface n.

[0172] In some embodiments, an inversion objective function is constructed based on a second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm. The objective function is constructed using the model parameter matrix and synthetic seismic records, specifically including:

[0173] The difference between the synthetic seismic record and the measured data is calculated to obtain the biased seismic record.

[0174] The difference between the model parameter matrix and the measured data is calculated to obtain the deviation model parameters.

[0175] An inversion objective function is constructed based on biased seismic records and biased model parameters.

[0176] The inversion objective function is expressed by the following formula:

[0177] J(m)=||(H+λ1I)Δm-J T d ref || 2 +λ2||Δm+m ref ||.

[0178] Where λ1 represents the damping coefficient, λ2 represents the sparse constraint coefficient, J represents the Jacobian matrix, and H represents the Hessian matrix H=J T J and I denote the identity matrix, d ref The variable represents the biased seismic record, and Δm represents the update amount of the model parameter matrix. ref This represents the parameters of the deviation model.

[0179] In some embodiments, the inversion objective function is solved and the model parameter matrix is ​​updated to obtain thin interlayer parameter information, specifically including:

[0180] Set the update amount of the model parameters, and iterate the inversion objective function based on the update amount of the model parameters.

[0181] In response to the inversion objective function approaching 0, the current model parameters are used as thin interlayer parameter information.

[0182] The update amount of the model parameters is expressed by the following formula:

[0183]

[0184] Where, d PP This indicates the PP wave synthesis record, d PS Indicates PS wave synthesis record, s PP This represents the measured data of the PP wave, s PS This represents the measured data of the PS wave.

[0185] This application also provides a thin interlayer inversion device, comprising:

[0186] The parameter construction module is used to construct the initial model parameter matrix based on preset seismic parameters.

[0187] The seismic calculation module calculates the reflection coefficient matrix of the subsurface medium based on the initial model parameter matrix, and uses the reflection coefficient matrix to perform forward modeling and synthesize seismic records.

[0188] The inversion module uses a second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt (LM) algorithm to construct the inversion objective function, which utilizes the model parameter matrix and synthetic seismic records.

[0189] The solver module is used to solve the inversion objective function and update the model parameter matrix to obtain thin interlayer parameter information.

[0190] In some embodiments, the seismic calculation module specifically includes:

[0191] The forward modeling module calculates the PP wave reflection coefficient and PS wave reflection coefficient matrix corresponding to the model parameters based on Kennett's second-order approximation formula.

[0192] The PP wave synthesis record and the PS wave synthesis record are obtained by calculation using the reflection coefficient matrix and the preset wavelet matrix.

[0193] Synthetic seismic records are calculated using the following formula:

[0194] d(m,θ,t)=W(t)R(m,θ,t).

[0195] Where d(m,θ,t) represents the synthetic seismic record, W(t) represents the wavelet matrix, R(m,θ,t) represents the reflection coefficient matrix, m represents the model parameter matrix, θ represents the reflection angle of the seismic wave, and t represents the time required for the seismic wave to travel from transmission to reception.

[0196] This application also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the method as described in any of the above.

[0197] This application also provides a non-transitory computer-readable storage medium that stores computer instructions for causing a computer to perform any of the methods described above.

[0198] Those skilled in the art should understand that the discussion of any of the above embodiments is merely exemplary and is not intended to imply that the scope of this application (including the claims) is limited to these examples; within the framework of this application, the technical features of the above embodiments or different embodiments can also be combined, the steps can be implemented in any order, and there are many other variations of different aspects of the embodiments of this application as described above, which are not provided in the details for the sake of brevity.

[0199] Additionally, to simplify the description and discussion, and to avoid obscuring the embodiments of this application, the well-known power / ground connections to integrated circuit (IC) chips and other components may or may not be shown in the provided drawings. Furthermore, the apparatus may be shown in block diagram form to avoid obscuring the embodiments of this application, and this also takes into account the fact that the details of the implementation of these block diagram apparatuses are highly dependent on the platform on which the embodiments of this application will be implemented (i.e., these details should be fully understood by those skilled in the art). While specific details (e.g., circuits) have been set forth to describe exemplary embodiments of this application, it will be apparent to those skilled in the art that the embodiments of this application can be implemented without these specific details or with variations thereof. Therefore, these descriptions should be considered illustrative rather than restrictive.

[0200] Although this application has been described in conjunction with specific embodiments thereof, many substitutions, modifications, and variations of these embodiments will be apparent to those skilled in the art from the foregoing description. For example, other memory architectures (e.g., dynamic RAM (DRAM)) may be used with the embodiments discussed.

[0201] The embodiments of this application are intended to cover all such substitutions, modifications, and variations that fall within the broad scope of the appended claims. Therefore, any omissions, modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the embodiments of this application should be included within the protection scope of this application.

Claims

1. A thin interlayer inversion method, comprising: An initial model parameter matrix is ​​constructed based on the parameterization of the underground medium; The reflection coefficient matrix of the subsurface medium is calculated based on the initial model parameter matrix, and forward modeling is performed using the reflection coefficient matrix to obtain a synthetic seismic record; wherein, the synthetic seismic record includes PP wave synthetic seismic record and PS wave synthetic seismic record; An inversion objective function is constructed based on the actual seismic records of the target subsurface medium and the synthetic seismic records; wherein, the actual seismic records include actual PP wave records and actual PS wave records; The inversion objective function is solved, and the initial parameter model is updated and iterated to obtain thin interlayer parameter information; The construction of the inversion objective function based on the actual seismic record and the synthetic seismic record of the target subsurface medium specifically includes: The difference between the synthetic seismic record and the measured data is calculated to obtain the biased seismic record; The difference between the model parameter matrix and the measured data is calculated to obtain the deviation model parameters; The inversion objective function is constructed based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt algorithm. The objective function is constructed using the model parameter matrix and the synthetic seismic record. Based on the Levenberg-Marquardt algorithm and with added sparsity constraints, the inversion objective function is expressed by the following formula: ； in, λ 1 represents the damping coefficient. λ 2 represents the sparsity constraint coefficient, J represents the Jacobian matrix, and H represents the Hessian matrix. I represents the identity matrix, d ref This indicates the aforementioned biased seismic record. This represents the update amount of the model parameter matrix. The parameters of the deviation model are represented; the J-expansion form under multi-wave joint inversion is: ； Solving the inversion objective function and updating the initial parameter model iteratively to obtain thin interlayer parameter information specifically includes: Find the solution when the objective function is sufficiently small to obtain the update amount of the model parameters, and then iterate over the initial model parameters; In response to the inversion objective function approaching 0, the current model parameters are used as the thin interlayer parameter information; The update amount of the model parameters is expressed by the following formula: ； Where, d PP This represents the PP wave synthesized seismic record, d PS This represents the PS wave synthesized seismic record, s PP This represents the measured data of the PP wave, s PS This represents the measured data of the PS wave.

2. The thin interlayer inversion method according to claim 1, wherein, The step of calculating the reflection coefficient matrix of the subsurface medium based on the initial model parameter matrix, and then using the reflection coefficient matrix to perform forward modeling and synthesize seismic records, specifically includes: The reflection coefficient matrix is ​​obtained by using the second-order approximation of the Kennett reflection coefficient formula based on the initial model parameter matrix. The PP wave composite seismic record and the PS wave composite seismic record are obtained by forward modeling based on the reflection coefficient matrix and the preset wavelet matrix. The synthetic seismic record is calculated using the following formula: ； in, This refers to the synthetic seismic record. Denotes the wavelet matrix, Represents the reflection coefficient matrix, This represents the model parameter matrix. Indicates the angle of reflection of seismic waves. This indicates the time required for a seismic wave to travel from transmission to reception.

3. The thin interlayer inversion method according to claim 2, wherein, The PP-wave composite seismic record and the PS-wave composite seismic record are obtained by forward modeling based on the reflection coefficient matrix and the preset wavelet matrix, specifically including: Calculate the PP wave reflection coefficient and PS wave reflection coefficient corresponding to the model parameters according to Kennett's second-order approximation formula; The PP wave composite seismic record and the PS wave composite seismic record are obtained by forward modeling based on the wave reflection coefficient and the preset wavelet matrix.

4. The thin interlayer inversion method according to claim 3, wherein, The calculation of the PP wave reflection coefficient and PS wave reflection coefficient corresponding to the model parameters according to Kennett's second-order approximation formula specifically includes: Construct a second-order approximate reflection coefficient matrix in the slowness-frequency domain based on Kennett's second-order approximation formula and the initial model parameter matrix; The second-order approximate reflection coefficient matrix is ​​calculated by the following formula: ； in, Represents angular frequency. p Indicates horizontal slowness. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n-1 for a downlink wave in the VTI medium. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for a downlink wave in the VTI medium. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the downlink wave. This represents the single-interface reflection coefficient matrix of the VTI medium at interface n for the upward wave. This represents the single-interface transmission coefficient matrix of the VTI medium at interface n for the upward wave. Let be the phase shift factor when the wave propagates from interface n-1 to interface n; Since the required reflection coefficient matrix needs to be in the angle-time domain, according to the inverse Fourier transform and Snell's law, the formula for calculating the total reflection coefficient matrix in the angle-time domain is: 。 5. A thin interlayer inversion device, comprising: The parameter construction module is used to construct the initial model parameter matrix based on the parameterization of the underground medium. The seismic calculation module is used to calculate the reflection coefficient matrix of the subsurface medium based on the initial model parameter matrix, and to perform forward modeling using the reflection coefficient matrix to obtain a synthetic seismic record; wherein, the synthetic seismic record includes a PP wave synthetic seismic record and a PS wave synthetic seismic record; The inversion module is used to construct an inversion objective function based on the actual seismic records of the target subsurface medium and the synthetic seismic records; wherein, the actual seismic records include actual PP wave records and actual PS wave records; The solution module is used to solve the inversion objective function, update and iterate the initial parameter model, and obtain thin interlayer parameter information; The construction of the inversion objective function based on the actual seismic record and the synthetic seismic record of the target subsurface medium specifically includes: The difference between the synthetic seismic record and the measured data is calculated to obtain the biased seismic record; The difference between the model parameter matrix and the measured data is calculated to obtain the deviation model parameters; The inversion objective function is constructed based on the second-order approximation of the Kennett reflection coefficient formula and the Levenberg-Marquardt algorithm. The objective function is constructed using the model parameter matrix and the synthetic seismic record. Based on the Levenberg-Marquardt algorithm and with added sparsity constraints, the inversion objective function is expressed by the following formula: ； in, λ 1 represents the damping coefficient. λ 2 represents the sparsity constraint coefficient, J represents the Jacobian matrix, and H represents the Hessian matrix. I represents the identity matrix, d ref This indicates the aforementioned biased seismic record. This represents the update amount of the model parameter matrix. The parameters of the deviation model are represented; the J-expansion form under multi-wave joint inversion is: ； Solving the inversion objective function and updating the initial parameter model iteratively to obtain thin interlayer parameter information specifically includes: Find the solution when the objective function is sufficiently small to obtain the update amount of the model parameters, and then iterate over the initial model parameters; In response to the inversion objective function approaching 0, the current model parameters are used as the thin interlayer parameter information; The update amount of the model parameters is expressed by the following formula: ； Where, d PP This represents the PP wave synthesized seismic record, d PS This represents the PS wave synthesized seismic record, s PP This represents the measured data of the PP wave, s PS This represents the measured data of the PS wave.

6. The thin interlayer inversion apparatus according to claim 5, wherein, The earthquake calculation module specifically includes: The forward modeling module is used to calculate the reflection coefficient matrix based on the initial model parameter matrix using a second-order approximation of the Kennett reflection coefficient formula. The PP wave composite seismic record and the PS wave composite seismic record are obtained by forward modeling based on the reflection coefficient matrix and the preset wavelet matrix. The synthetic seismic record is calculated using the following formula: ； in, This refers to the synthetic seismic record. Denotes the wavelet matrix, Represents the reflection coefficient matrix, This represents the model parameter matrix. Indicates the angle of reflection of seismic waves. This indicates the time required for a seismic wave to travel from transmission to reception.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the method as claimed in any one of claims 1 to 4.

8. A non-transitory computer-readable storage medium storing computer instructions for causing a computer to perform the method of any one of claims 1 to 4.

Images