Frequency-dependent fluid mobility inversion method and device

Abstract

The invention discloses a frequency-dependent fluid mobility inversion method and device, belongs to the technical field of seismic exploration, and solves the problem that a fluid mobility result extracted by an existing time-frequency analysis method is difficult to directly reflect real physical distribution of a reservoir. The method comprises the following steps: acquiring seismic data including multiple seismic post-stack records and relative wave impedance; the number of frequency bands needing to participate in inversion and the frequency range corresponding to each frequency band are set, continuous wavelet transform is carried out on each seismic post-stack record, and a continuous wavelet transform result is inversely transformed to a time domain to obtain seismic frequency division data needing to participate in inversion; calculating a fluid mobility coefficient of the reservoir according to a logging curve obtained by logging, and interpolating the fluid mobility coefficient into a low-frequency model; and based on the relative wave impedance and the L2 regularization constraint, establishing a fluid mobility inversion equation, substituting the seismic frequency division data, the low-frequency model and the relative wave impedance into the established fluid mobility inversion equation, and solving a mobility inversion result recorded after each seismic stack by using gradient descent.

Description

Technical Field

[0001] A method and device for frequency-dependent fluid mobility inversion, which are used to obtain frequency-divided data by continuous wavelet transform and invert fluid mobility by using L2 norm regularization, belong to the field of seismic exploration technology. Background Art

[0002] Fluid mobility is a key geophysical parameter used in seismic exploration to characterize the flow capacity of subsurface fluids (such as oil, gas, and water) within a reservoir. Defined as the ratio of formation permeability to fluid viscosity, it effectively reflects the pore structure of the reservoir rock skeleton and the fluid's mobility. In oil and gas exploration, fluid mobility can be used to identify oil and gas reservoirs and assist in reservoir prediction, and is widely used in the exploration and development of both conventional and unconventional oil and gas resources.

[0003] Existing technologies primarily rely on time-frequency analysis methods to extract mobility attributes, such as short-time Fourier transform (STFT), continuous wavelet transform (CWT), generalized S-transform (GST), and empirical mode decomposition (EMD). These methods indirectly calculate mobility attributes by analyzing the characteristics of seismic signals at different frequency and time scales. However, existing methods suffer from the following technical issues during application:

[0004] 1. The mobility extraction method of time-frequency analysis is more biased towards the interface characteristics of seismic signals and cannot directly reflect the actual physical distribution of the reservoir;

[0005] 2. Strong dependence on seismic data quality: Time-frequency analysis methods are sensitive to parameters such as the signal-to-noise ratio and resolution of seismic data. Fluctuations in data quality may lead to unstable mobility calculation results.

[0006] 3. Complex parameter adjustment: Different time-frequency analysis methods involve multiple parameters (such as window length, decomposition order, mother wavelet type, etc.), which need to be adjusted according to different geological conditions. The operation is cumbersome and the computational cost is high. Summary of the Invention

[0007] In response to the above research problems, the purpose of the present invention is to provide a method and device for frequency-dependent fluid mobility inversion to solve the problem that the mobility extraction of existing time-frequency analysis methods is more biased towards the interface characteristics of seismic signals and is difficult to directly reflect the true physical distribution of the reservoir.

[0008] In order to achieve the above object, the present invention adopts the following technical solutions:

[0009] A method for frequency-dependent fluid mobility inversion comprises the following steps:

[0010] Step 1: Obtain seismic data including multi-channel seismic post-stack records and relative wave impedance;

[0011] Step 2: Set the number of frequency bands required for inversion and the frequency range corresponding to each frequency band, and perform continuous wavelet transform on each seismic post-stack record. According to the set number of frequency bands and frequency band range, the continuous wavelet transform result is inversely transformed to the time domain to obtain the seismic frequency data required for inversion;

[0012] Step 3: Calculate the fluid mobility coefficient of the reservoir based on the logging curves obtained from well logging and interpolate it into a low-frequency model;

[0013] Step 4: Based on the relative wave impedance and L2 regularization constraints, establish the fluid mobility inversion equation. Substitute the seismic frequency-divided data, low-frequency model, and relative wave impedance into the established fluid mobility inversion equation, and use gradient descent to solve the mobility inversion results of each seismic post-stack record.

[0014] Furthermore, the specific steps of step 2 are:

[0015] Step 2.1: Set the number of frequency bands to be inverted, and set the frequency range of each frequency band, including the start frequency and end frequency;

[0016] Step 2.2: Perform continuous wavelet transform on each post-stack earthquake record x(t) corresponding to the frequency band to obtain the continuous wavelet transform result of each post-stack earthquake record x(t). The formula of continuous wavelet transform is:

[0017]

[0018]

[0019] Where W(a, b) represents the continuous wavelet transform result of each post-stack earthquake recording signal x(t), a represents the wavelet scale factor, b represents the translation factor, which is used to control the time position, and ψ a，b (t) represents the mother wavelet after translation by b and scaling by a at time t, * represents the complex conjugate, and ψ represents the mother wavelet;

[0020] Step 2.3: Transform the continuous wavelet transform results of the required frequency band into the time domain to obtain the seismic frequency data required for inversion. The formula is:

[0021]

[0022] The corresponding relationship between frequency and scale after continuous wavelet transform is expressed as:

[0023]

[0024]

[0025] W filtered(a, b) = W(a, b)·E(a) (6)

[0026]

[0027] Where s(t) is the inverted seismic frequency data, W filtered (a, b) means in a min <a<a max The new wavelet coefficient result constructed by the corresponding frequency band, C ψ represents the volume constant of the wavelet, f represents the frequency, F c represents the center frequency of the mother wavelet, [f min , f max ] represents the frequency band range, f min <f<f max , a max Indicates the end point of the frequency band f max The corresponding scale, a min Indicates the starting point of the frequency band f min The corresponding scale, E(a) represents the Dirac function.

[0028] Furthermore, the specific steps of step 3 are:

[0029] Step 3.1: Obtain the well logging curve based on the well logging curve, based on the P-wave velocity of the reservoir at the i-th sampling point , shear wave velocity Density ρ i and porosity φ i Calculate the fluid mobility coefficient of each sampling point using the following formula:

[0030]

[0031]

[0032]

[0033]

[0034]

[0035]

[0036]

[0037]

[0038]

[0039]

[0040]

[0041]

[0042] Where C(i) represents the fluid mobility coefficient of the i-th sampling point in the reservoir, R1(i) represents the first-order term of the reflection coefficient of the i-th sampling point in the reservoir, represents the fluid density at the i-th sampling point in the reservoir, Z i represents the impedance of the i-th sampling point in the reservoir, represents the equivalent P-wave velocity of the i-th sampling point in the reservoir, in dry rock M i is the bulk modulus or longitudinal wave modulus of the i-th sampling point in the reservoir, represents the bulk modulus of dry rock at the i-th sampling point in the reservoir, μ i represents the second Lame coefficient of the i-th sampling point in the reservoir, and represents the dimensionless coefficient of the i-th sampling point in the reservoir, represents the propagation velocity of the wave in the fluid at the i-th sampling point in the reservoir, represents the compressibility of the fluid at the i-th sampling point in the reservoir, and denote the elastic moduli of the skeleton and fluid at the i-th sampling point in the reservoir, represents the bulk modulus of the rock matrix at the i-th sampling point in the reservoir, which is obtained using the LRM linear fitting method. i 、D i 、P i , Q i It is an intermediate variable used to simplify the formula during the calculation process;

[0043] Step 3.2: Interpolate the fluid mobility coefficient of each sampling point in the reservoir into a low-frequency model.

[0044] Furthermore, the specific steps of step 4 are:

[0045] Step 4.1: Establish the fluid mobility inversion equation based on relative wave impedance and L2 regularization constraint. The specific steps are as follows:

[0046] According to Silin, when a fast longitudinal wave is incident vertically, the reflection coefficient at the interface of a viscous fluid-porous medium can be described as follows in the low-frequency asymptotic form:

[0047]

[0048] Where R is the longitudinal wave reflection coefficient, ε is a small dimensionless parameter, |ε|=|jλω|, j is an imaginary unit, ω is the frequency, ρf is the density of viscous fluid, k is the permeability of porous rock, η is the viscosity coefficient of fluid, Z1 and Z2 are the longitudinal wave impedances of the upper and lower media respectively, λ is the intermediate variable, and R1 represents the first-order term of the reflection coefficient;

[0049] Write equation (20) as a function of time and frequency, and take the real part to get:

[0050]

[0051] in, is the mobility coefficient at time t, which means that in the inversion process, the low-frequency model established by calculating the fluid mobility coefficient of all sampling points in the logging reservoir is replaced by K(t). It represents the fluid mobility value of all sampling points of the logging reservoir at time t;

[0052] Write (21) in matrix form and expand it according to frequency to obtain the following equations:

[0053]

[0054] Among them, I n×n is the n-dimensional identity matrix, ω m represents the mth frequency sampling point, is an n-dimensional diagonal matrix;

[0055] Multiply both ends of equation (22) by the wavelet matrix to obtain:

[0056]

[0057] Formula (23) can also be simplified into the following form:

[0058] G·X0=d (24)

[0059] Where W is the wavelet matrix, d(t, ω m ) is the earthquake frequency data according to frequency, G represents X0 represents d represents

[0060] If the seismic wave is incident at a small angle, the relationship between the reflection coefficient and the longitudinal wave impedance is in the form of time and frequency:

[0061]

[0062] Where: I P (t, f) is the longitudinal wave impedance in the form of time t and frequency f, I P (t0, f) is the initial impedance value,

[0063] Combining equations (22) and (25) and expanding them according to frequency, we get the model constraint equation:

[0064]

[0065] in, It is expanded into frequency-dependent impedance, D is the summation matrix, and Equation (26) is simplified into the following form:

[0066] C′·X0=ε (27)

[0067] In the formula, C′ represents X0 represents ε represents

[0068] Combining equations (27) and (24), we construct the inversion objective function:

[0069]

[0070] Introducing the L2 norm into the inversion objective function yields a new objective function:

[0071]

[0072] Where α is the impedance model Parameters, β is the regularization coefficient;

[0073] Step 4.2: Gradient descent method is used to solve the fluid mobility, that is, to solve equation (29) based on the gradient descent method. The specific steps are as follows:

[0074] Based on the new objective function (29), the gradient formula is:

[0075]

[0076] Where G T represents the transpose of the G matrix, C′ T It is represented as the transpose of the C′ matrix, and E is represented as the identity matrix;

[0077] The gradient of the first round is obtained based on the gradient formula: If satisfied Then the mobility inversion results of each earthquake post-stack record are obtained Otherwise, the update using gradient descent is Make another judgment.

[0078]

[0079] Where η is the learning rate; this formula indicates that when the second norm of the gradient is less than the threshold ε, the gradient descent algorithm stops, indicating that it is close to the optimal solution.

[0080] A device for frequency-dependent fluid mobility inversion, comprising:

[0081] Acquisition module: acquires seismic data including multi-channel seismic post-stack records and relative wave impedance;

[0082] Seismic frequency data acquisition module: Set the number of frequency bands required for inversion and the frequency range corresponding to each frequency band, and perform continuous wavelet transform on each seismic post-stack record. According to the set number of frequency bands and frequency band range, the continuous wavelet transform result is inversely transformed into the time domain to obtain the seismic frequency data required for inversion;

[0083] Interpolation module: Calculates the fluid mobility coefficient of the reservoir based on the logging curves obtained from well logging and interpolates it into a low-frequency model;

[0084] Inversion module: Based on relative wave impedance and L2 regularization constraints, the fluid mobility inversion equation is established. The seismic frequency-dividing data, low-frequency model and relative wave impedance are introduced into the established fluid mobility inversion equation, and the mobility inversion results of each seismic post-stack record are solved using gradient descent.

[0085] Furthermore, the specific implementation steps of the earthquake frequency division data acquisition module are as follows:

[0086] Step 2.1: Set the number of frequency bands to be inverted, and set the frequency range of each frequency band, including the start frequency and end frequency;

[0087] Step 2.2: Perform continuous wavelet transform on each post-stack earthquake record x(t) corresponding to the frequency band to obtain the continuous wavelet transform result of each post-stack earthquake record x(t). The formula of continuous wavelet transform is:

[0088]

[0089]

[0090] Where W(a, b) represents the continuous wavelet transform result of each post-stack earthquake recording signal x(t), a represents the wavelet scale factor, b represents the translation factor, which is used to control the time position, and ψ a，b (t) represents the mother wavelet after translation by b and scaling by a at time t, * represents the complex conjugate, and ψ represents the mother wavelet;

[0091] Step 2.3: Transform the continuous wavelet transform results of the required frequency band into the time domain to obtain the seismic frequency data required for inversion. The formula is:

[0092]

[0093] The corresponding relationship between frequency and scale after continuous wavelet transform is expressed as:

[0094]

[0095]

[0096] W filtered (a, b) = W(a, b)·E(a) (6)

[0097]

[0098] Where s(t) is the inverted seismic frequency data, W filtered (a, b) means in a min <a<a max The new wavelet coefficient result constructed by the corresponding frequency band, C ψ represents the volume constant of the wavelet, f represents the frequency, F c represents the center frequency of the mother wavelet, [f min , f max ] represents the frequency band range, f min <f<f max , a max Indicates the end point of the frequency band f max The corresponding scale, a min Indicates the starting point of the frequency band f min The corresponding scale, E(a) represents the Dirac function.

[0099] Furthermore, the specific implementation steps of the interpolation module are:

[0100] Step 3.1: Obtain the well logging curve based on the well logging curve, based on the P-wave velocity of the reservoir at the i-th sampling point , shear wave velocity Density ρ i and porosity φ i Calculate the fluid mobility coefficient of each sampling point using the following formula:

[0101]

[0102]

[0103]

[0104]

[0105]

[0106]

[0107]

[0108]

[0109]

[0110]

[0111]

[0112]

[0113] Where C(i) represents the fluid mobility coefficient of the i-th sampling point in the reservoir, R1(i) represents the first-order term of the reflection coefficient of the i-th sampling point in the reservoir, represents the fluid density at the i-th sampling point in the reservoir, Z i represents the impedance of the i-th sampling point in the reservoir, represents the equivalent P-wave velocity of the i-th sampling point in the reservoir, in dry rock M i is the bulk modulus or longitudinal wave modulus of the i-th sampling point in the reservoir, represents the bulk modulus of dry rock at the i-th sampling point in the reservoir, μ i represents the second Lame coefficient of the i-th sampling point in the reservoir, and represents the dimensionless coefficient of the i-th sampling point in the reservoir, represents the propagation velocity of the wave in the fluid at the i-th sampling point in the reservoir, represents the compressibility of the fluid at the i-th sampling point in the reservoir, and denote the elastic moduli of the skeleton and fluid at the i-th sampling point in the reservoir, represents the bulk modulus of the rock matrix at the i-th sampling point in the reservoir, which is obtained using the LRM linear fitting method. i 、D i 、P i , Q i It is an intermediate variable used to simplify the formula during the calculation process;

[0114] Step 3.2: Interpolate the fluid mobility coefficient of each sampling point in the reservoir into a low-frequency model.

[0115] Furthermore, the specific implementation steps of the inversion module are:

[0116] Step 4.1: Establish the fluid mobility inversion equation based on relative wave impedance and L2 regularization constraint. The specific steps are as follows:

[0117] According to Silin, when a fast longitudinal wave is incident vertically, the reflection coefficient at the interface of a viscous fluid-porous medium can be described as follows in the low-frequency asymptotic form:

[0118]

[0119] Where R is the longitudinal wave reflection coefficient, ε is a small dimensionless parameter, |ε|=|jλω|, j is an imaginary unit, ω is the frequency, ρ f is the density of viscous fluid, k is the permeability of porous rock, η is the viscosity coefficient of fluid, Z1 and Z2 are the longitudinal wave impedances of the upper and lower media respectively, λ is the intermediate variable, and R1 represents the first-order term of the reflection coefficient;

[0120] Write equation (20) as a function of time and frequency, and take the real part to get:

[0121]

[0122] in, is the mobility coefficient at time t, which means that in the inversion process, the low-frequency model established by calculating the fluid mobility coefficient of all sampling points in the logging reservoir is replaced by K(t). It represents the fluid mobility value of all sampling points of the logging reservoir at time t;

[0123] Write (21) in matrix form and expand it according to frequency to obtain the following equations:

[0124]

[0125] Among them, I n×n is the n-dimensional identity matrix, ω m represents the mth frequency sampling point, is an n-dimensional diagonal matrix;

[0126] Multiply both ends of equation (22) by the wavelet matrix to obtain:

[0127]

[0128] Formula (23) can also be simplified into the following form:

[0129] G·X0=d (24)

[0130] Where W is the wavelet matrix, d(t, ω m ) is the earthquake frequency data according to frequency, G represents X0 represents d represents

[0131] If the seismic wave is incident at a small angle, the relationship between the reflection coefficient and the longitudinal wave impedance is in the form of time and frequency:

[0132]

[0133] Where: I P (t, f) is the longitudinal wave impedance in the form of time t and frequency f, I P (t0, f) is the initial impedance value,

[0134] Combining equations (22) and (25) and expanding them according to frequency, we get the model constraint equation:

[0135]

[0136] in, It is expanded into frequency-dependent impedance, D is the summation matrix, and Equation (26) is simplified into the following form:

[0137] C′·X0=ε (27)

[0138] Where C′ represents X0 represents ε represents

[0139] Combining equations (27) and (24), we construct the inversion objective function:

[0140]

[0141] Introducing the L2 norm into the inversion objective function yields a new objective function:

[0142]

[0143] Where α is the impedance model Parameters, β is the regularization coefficient;

[0144] Step 4.2: Gradient descent method is used to solve the fluid mobility, that is, to solve equation (29) based on the gradient descent method. The specific steps are as follows:

[0145] Based on the new objective function (29), the gradient formula is:

[0146]

[0147] Where G T represents the transpose of the G matrix, C′ T It is represented as the transpose of the C′ matrix, and E is represented as the identity matrix;

[0148] The gradient of the first round is obtained based on the gradient formula: If satisfied Then the mobility inversion results of each earthquake post-stack record are obtained Otherwise, the update using gradient descent is Make another judgment.

[0149]

[0150] Where η is the learning rate; this formula indicates that when the second norm of the gradient is less than the threshold ε, the gradient descent algorithm stops, indicating that it is close to the optimal solution.

[0151] Compared with the prior art, the present invention has the following beneficial effects:

[0152] First, the present invention uses a seismic inversion model to directly model mobility attributes as unknown variables and obtains them through optimization, effectively bypassing the limitations of traditional methods that are based solely on waveform interfaces. The extracted results are more physically consistent and can more realistically reflect fluid changes within the reservoir. This avoids the problem that mobility extraction using time-frequency analysis methods is more biased towards the interface characteristics of seismic signals and cannot directly reflect the true physical distribution of the reservoir.

[0153] Second, the present invention introduces an L2 norm regularized mobility inversion method to improve the stability and robustness of the solution. This avoids the sensitivity of the time-frequency analysis method to parameters such as the signal-to-noise ratio and resolution of seismic data, which may lead to unstable mobility calculation results due to fluctuations in data quality.

[0154] 3. The present invention only introduces a single regularization parameter, which greatly simplifies the parameter setting process, making it easy to operate and low in computational cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0155] Figure 1 It is a schematic diagram of the flow chart of the present invention.

[0156] Figure 2 Schematic diagram of the impedance model (ie, relative wave impedance) in the present invention.

[0157] Figure 3 Schematic diagram of multi-channel seismic post-stack recording in the present invention.

[0158] Figure 4 Schematic diagram of the mobility extraction results calculated by traditional methods.

[0159] Figure 5 Schematic diagram of the mobility inversion result in the present invention. DETAILED DESCRIPTION

[0160] The present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0161] The following combination Figure 1-5 The present invention is described in detail.

[0162] By utilizing a post-stack mobility inversion method, the present invention effectively overcomes the problem of mobility extraction methods being biased towards seismic signal interface characteristics, thereby more accurately reflecting the true physical distribution of the reservoir. Seismic data comprising multiple seismic post-stack records and relative wave impedance are obtained. The number of frequency bands to be inverted and the frequency range corresponding to each frequency band are set. Based on this, each post-stack seismic record is extracted and a continuous wavelet transform is performed on each post-stack seismic record. Based on the set number of frequency bands and frequency range, the continuous wavelet transform is inversely transformed to the time domain to obtain the seismic frequency-dividing data required for inversion. The mobility coefficient is calculated based on the well logging curve and interpolated into a low-frequency model. Based on relative wave impedance and L2 regularization constraints, a fluid mobility inversion equation is established. The seismic frequency-dividing data, low-frequency model, and relative wave impedance are introduced into the fluid mobility inversion equation, and the mobility inversion result for each channel is solved using gradient descent. The mobility inversion method can penetrate deep into the reservoir, providing a mobility distribution that is more consistent with geological reality, greatly improving the prediction accuracy and reliability of oil and gas reservoirs.

[0163] Specifically, the present invention provides a method for frequency-dependent fluid mobility inversion, such as Figure 1 As shown, the following steps are included:

[0164] Step 1: Obtain seismic data including multi-channel seismic post-stack records and relative wave impedance, such as Figure 3 Shown is a multi-channel seismic post-stack record;

[0165] Step 2: Set the number of frequency bands required for inversion and the frequency range corresponding to each frequency band, and perform continuous wavelet transform on each seismic post-stack record. According to the set number of frequency bands and frequency band range, the continuous wavelet transform result is inversely transformed to the time domain to obtain the seismic frequency data required for inversion;

[0166] The specific steps are:

[0167] Step 2.1: Set the number of frequency bands to be inverted, and set the frequency range of each frequency band, including the start frequency and end frequency;

[0168] Step 2.2: Perform continuous wavelet transform on each post-stack earthquake record x(t) corresponding to the frequency band to obtain the continuous wavelet transform result of each post-stack earthquake record x(t). The formula of continuous wavelet transform is:

[0169]

[0170]

[0171] Where W(a, b) represents the continuous wavelet transform result of each post-stack earthquake recording signal x(t), a represents the wavelet scale factor, b represents the translation factor, which is used to control the time position, and ψa，b (t) represents the mother wavelet after translation by b and scaling by a at time t, * represents the complex conjugate, and ψ represents the mother wavelet;

[0172] Step 2.3: Transform the continuous wavelet transform results of the required frequency band into the time domain to obtain the seismic frequency data required for inversion. The formula is:

[0173]

[0174] The corresponding relationship between frequency and scale after continuous wavelet transform is expressed as:

[0175]

[0176]

[0177] W filtered (a, b) = W(a, b)·E(a) (6)

[0178]

[0179] Where s(t) is the inverted seismic frequency data, W filtered (a, b) means in a min <a<a max The new wavelet coefficient result constructed by the corresponding frequency band, C ψ represents the volume constant of the wavelet, f represents the frequency, F c represents the center frequency of the mother wavelet, [f min , f max ] represents the frequency band range, f min <f<f max , a max Indicates the end point of the frequency band f max The corresponding scale, a min Indicates the starting point of the frequency band f min The corresponding scale, E(a) represents the Dirac function.

[0180] Step 3: Calculate the fluid mobility coefficient of the reservoir based on the logging curves obtained from well logging and interpolate it into a low-frequency model;

[0181] The specific steps are:

[0182] Step 3.1: Obtain the well logging curve based on the well logging curve, based on the P-wave velocity of the reservoir at the i-th sampling point , shear wave velocity Density ρ i and porosity φ i Calculate the fluid mobility coefficient of each sampling point using the following formula:

[0183]

[0184]

[0185]

[0186]

[0187]

[0188]

[0189]

[0190]

[0191]

[0192] Where C(i) represents the fluid mobility coefficient of the i-th sampling point in the reservoir, R1(i) represents the first-order term of the reflection coefficient of the i-th sampling point in the reservoir, represents the fluid density at the i-th sampling point in the reservoir, Z i represents the impedance of the i-th sampling point in the reservoir, represents the equivalent P-wave velocity of the i-th sampling point in the reservoir, in dry rock M i is the bulk modulus or longitudinal wave modulus of the i-th sampling point in the reservoir, represents the bulk modulus of dry rock at the i-th sampling point in the reservoir, μ i represents the second Lame coefficient of the i-th sampling point in the reservoir, and represents the dimensionless coefficient of the i-th sampling point in the reservoir, represents the propagation velocity of the wave in the fluid at the i-th sampling point in the reservoir, represents the compressibility of the fluid at the i-th sampling point in the reservoir, and denote the elastic moduli of the skeleton and fluid at the i-th sampling point in the reservoir, represents the bulk modulus of the rock matrix at the i-th sampling point in the reservoir, which is obtained using the LRM linear fitting method. i 、D i 、P i , Q i It is an intermediate variable used to simplify the formula during the calculation process;

[0193] Step 3.2: Interpolate the fluid mobility coefficient of each sampling point in the reservoir into a low-frequency model.

[0194] Step 4: Based on the relative wave impedance and L2 regularization constraints, establish the fluid mobility inversion equation, bring the seismic frequency data, low-frequency model and relative wave impedance into the established fluid mobility inversion equation, and use gradient descent to solve the mobility inversion results of each seismic post-stack record, such as Figure 5 shown.

[0195] Step 4.1: Establish the fluid mobility inversion equation based on relative wave impedance and L2 regularization constraint. The specific steps are as follows:

[0196] According to Silin, when a fast longitudinal wave is incident vertically, the reflection coefficient at the interface of a viscous fluid-porous medium can be described as follows in the low-frequency asymptotic form:

[0197]

[0198] Where R is the longitudinal wave reflection coefficient, ε is a small dimensionless parameter, |ε|=|jλω|, j is an imaginary unit, ω is the frequency, ρ f is the density of viscous fluid, k is the permeability of porous rock, η is the viscosity coefficient of fluid, Z1 and Z2 are the longitudinal wave impedances of the upper and lower media respectively, λ is the intermediate variable, and R1 represents the first-order term of the reflection coefficient;

[0199] Write equation (20) as a function of time and frequency, and take the real part to get:

[0200]

[0201] in, is the mobility coefficient at time t, which means that in the inversion process, the low-frequency model established by calculating the fluid mobility coefficient of all sampling points in the logging reservoir is replaced by K(t). It represents the fluid mobility value of all sampling points of the logging reservoir at time t;

[0202] Write (21) in matrix form and expand it according to frequency to obtain the following equations:

[0203]

[0204] Among them, I n×n is the n-dimensional identity matrix, ω m represents the mth frequency sampling point, is an n-dimensional diagonal matrix;

[0205] Multiply both ends of equation (22) by the wavelet matrix to obtain:

[0206]

[0207] Formula (23) can also be simplified into the following form:

[0208] G·X0=d (24)

[0209] Where W is the wavelet matrix, d(t, ω m ) is the earthquake frequency data according to frequency, G represents X0 represents d represents

[0210] If the seismic wave is incident at a small angle, the relationship between the reflection coefficient and the longitudinal wave impedance is in the form of time and frequency:

[0211]

[0212] Where: I P (t, f) is the longitudinal wave impedance in the form of time t and frequency f, I P (t0, f) is the initial impedance value,

[0213] Combining equations (22) and (25) and expanding them according to frequency, we get the model constraint equation:

[0214]

[0215] in, It is expanded into frequency-dependent impedance, D is the summation matrix, and Equation (26) is simplified into the following form:

[0216] C′·X0=ε (27)

[0217] Where C′ represents X0 represents ε represents

[0218] Combining equations (27) and (24), we construct the inversion objective function:

[0219]

[0220] Introducing the L2 norm into the inversion objective function yields a new objective function:

[0221]

[0222] Where α is the impedance model Parameters such as Figure 2 As shown, β is the regularization coefficient;

[0223] Step 4.2: Gradient descent method is used to solve the fluid mobility, that is, to solve equation (29) based on the gradient descent method. The specific steps are as follows:

[0224] Based on the new objective function (29), the gradient formula is:

[0225]

[0226] Where G T represents the transpose of the G matrix, C′ T It is represented as the transpose of the C′ matrix, and E is represented as the identity matrix;

[0227] The gradient of the first round is obtained based on the gradient formula: If satisfied Then the mobility inversion results of each earthquake post-stack record are obtained Otherwise, the update using gradient descent is Make another judgment.

[0228]

[0229] Where η is the learning rate; this formula indicates that when the second norm of the gradient is less than the threshold ε, the gradient descent algorithm stops, indicating that it is close to the optimal solution.

[0230] Effect analysis: Figure 2-5 As shown, Figure 4 The mobility results extracted using the time-frequency analysis method show that they are more inclined to the interface characteristics of the seismic signal and are difficult to directly reflect the real physical distribution of the reservoir. Figure 5 The fluid mobility inversion results obtained in this study reflect fluid locations that are closer to the geological model. The mobility inversion method can penetrate deep into the reservoir, providing a mobility distribution that is more consistent with geological reality, proving the correctness of this method. The fluid mobility obtained using this method has higher resolution than traditional methods, making it potentially applicable to oil and gas reservoir identification. This improves the recognition accuracy of fluid detection and can more accurately indicate the location of the reservoir.

[0231] The above are only representative embodiments of the present invention in many specific application scopes and do not constitute any limitation on the protection scope of the present invention. Any technical solutions formed by transformation or equivalent replacement fall within the scope of protection of the present invention.

Claims

1. A method for frequency-dependent fluid mobility inversion, characterized in that: The steps include: Step 1: Obtain seismic data including multi-channel seismic post-stack records and relative wave impedance; Step 2: Set the number of frequency bands required for inversion and the frequency range corresponding to each frequency band, and perform continuous wavelet transform on each seismic post-stack record. According to the set number of frequency bands and frequency band range, the continuous wavelet transform result is inversely transformed to the time domain to obtain the seismic frequency data required for inversion; Step 3: Calculate the fluid mobility coefficient of the reservoir based on the logging curves obtained from well logging and interpolate it into a low-frequency model; Step 4: Based on the relative wave impedance and L2 regularization constraints, establish the fluid mobility inversion equation. Substitute the seismic frequency-divided data, low-frequency model, and relative wave impedance into the established fluid mobility inversion equation, and use gradient descent to solve the mobility inversion results of each seismic post-stack record.

2. A method for frequency-dependent fluid mobility inversion according to claim 1, characterized in that: The specific steps of step 2 are: Step 2.1: Set the number of frequency bands to be inverted, and set the frequency range of each frequency band, including the start frequency and end frequency; Step 2.2: Perform continuous wavelet transform on each post-stack earthquake record x(t) corresponding to the frequency band to obtain the continuous wavelet transform result of each post-stack earthquake record x(t). The formula of continuous wavelet transform is: Where W(a, b) represents the continuous wavelet transform result of each post-stack seismic recording signal x(t), a represents the wavelet scale factor, b represents the translation factor, which is used to control the time position, and ψ a，b (t) represents the mother wavelet after translation by b and scaling by a at time t, * represents the complex conjugate, and ψ represents the mother wavelet; Step 2.3: Transform the continuous wavelet transform results of the required frequency band into the time domain to obtain the seismic frequency data required for inversion. The formula is: The corresponding relationship between frequency and scale after continuous wavelet transform is expressed as: W filtered (a,b)=W(a,b)·E(a) (6) Where s(t) is the inverted seismic frequency data, W filtered (a, b) means in a min <a max The new wavelet coefficient result constructed by the corresponding frequency band, C ψ represents the volume constant of the wavelet, f represents the frequency, F c represents the center frequency of the mother wavelet, [f min , f max ] represents the frequency band range, f min <f<f max , a max Indicates the end point of the frequency band f max The corresponding scale, a min Indicates the starting point of the frequency band f min The corresponding scale, E(a) represents the Dirac function.​ 3. The method for frequency-dependent fluid mobility inversion according to claim 2, characterized in that: The specific steps of step 3 are: Step 3.1: Obtain the well logging curve based on the well logging curve, based on the P-wave velocity of the reservoir at the i-th sampling point , shear wave velocity Density ρ i and porosity φ i Calculate the fluid mobility coefficient of each sampling point using the following formula: Where C(i) represents the fluid mobility coefficient of the i-th sampling point in the reservoir, R1(i) represents the first-order term of the reflection coefficient of the i-th sampling point in the reservoir, represents the fluid density at the i-th sampling point in the reservoir, Z i represents the impedance of the i-th sampling point in the reservoir, represents the equivalent P-wave velocity of the i-th sampling point in the reservoir, in dry rock M i is the bulk modulus or longitudinal wave modulus of the i-th sampling point in the reservoir, represents the bulk modulus of dry rock at the i-th sampling point in the reservoir, μ i represents the second Lame coefficient of the i-th sampling point in the reservoir, and represents the dimensionless coefficient of the i-th sampling point in the reservoir, represents the propagation velocity of the wave in the fluid at the i-th sampling point in the reservoir, represents the compressibility of the fluid at the i-th sampling point in the reservoir, and denote the elastic moduli of the skeleton and fluid at the i-th sampling point in the reservoir, represents the bulk modulus of the rock matrix at the i-th sampling point in the reservoir, which is obtained using the LRM linear fitting method. i 、D i 、P i , Q i It is an intermediate variable used to simplify the formula during the calculation process; Step 3.2: Interpolate the fluid mobility coefficient of each sampling point in the reservoir into a low-frequency model.

4. The method for frequency-dependent fluid mobility inversion according to claim 3, characterized in that: The specific steps of step 4 are: Step 4.1: Establish the fluid mobility inversion equation based on relative wave impedance and L2 regularization constraint. The specific steps are as follows: According to Silin, when a fast longitudinal wave is incident vertically, the reflection coefficient at the interface of a viscous fluid-porous medium can be described as follows in the low-frequency asymptotic form: Where R is the longitudinal wave reflection coefficient, ε is a small dimensionless parameter, |ε| = |jλω|, j is an imaginary unit, ω is the frequency, ρ f is the density of viscous fluid, k is the permeability of porous rock, η is the viscosity coefficient of fluid, Z1 and Z2 are the longitudinal wave impedances of the upper and lower media respectively, λ is the intermediate variable, and R1 represents the first-order term of the reflection coefficient; Write equation (20) as a function of time and frequency, and take the real part to get: in, is the mobility coefficient at time t, which means that in the inversion process, the low-frequency model established by calculating the fluid mobility coefficient of all sampling points in the logging reservoir is replaced by K(t). It represents the fluid mobility value of all sampling points of the logging reservoir at time t; Write (21) in matrix form and expand it according to frequency to obtain the following equations: Among them, I n×n is the n-dimensional identity matrix, ω m represents the mth frequency sampling point, is an n-dimensional diagonal matrix; Multiply both ends of equation (22) by the wavelet matrix to obtain: Formula (23) can also be simplified into the following form: G·X0=d (24) Where W is the wavelet matrix, d(t, ω m ) is the earthquake frequency data according to frequency, G represents X0 represents d represents If the seismic wave is incident at a small angle, the relationship between the reflection coefficient and the longitudinal wave impedance is in the form of time and frequency: Where: I P (t, f) is the longitudinal wave impedance in the form of time t and frequency f, I P (t0, f) is the initial impedance value, Combining equations (22) and (25) and expanding them according to frequency, we get the model constraint equation: in, It is expanded into frequency-dependent impedance, D is the summation matrix, and Equation (26) is simplified into the following form: C·X0=ε (27) Where C represents X0 represents ε represents Combining equations (27) and (24), we construct the inversion objective function: Introducing the L2 norm into the inversion objective function yields a new objective function: Where α is the impedance model Parameters, β is the regularization coefficient; Step 4.2: Gradient descent method is used to solve the fluid mobility, that is, to solve equation (29) based on the gradient descent method. The specific steps are as follows: Based on the new objective function (29), the gradient formula is: Where G T represents the transpose of the G matrix, C′ T It is represented as the transpose of the C′ matrix, and E is represented as the identity matrix; The gradient of the first round is obtained based on the gradient formula: If satisfied Then the mobility inversion results of each earthquake post-stack record are obtained Otherwise, the update using gradient descent is Then make a judgment; Where η is the learning rate; this formula indicates that when the second norm of the gradient is less than the threshold ε, the gradient descent algorithm stops, indicating that it is close to the optimal solution.

5. A device for frequency-dependent fluid mobility inversion, characterized in that: include: Acquisition module: acquires seismic data including multi-channel seismic post-stack records and relative wave impedance; Seismic frequency data acquisition module: Set the number of frequency bands required for inversion and the frequency range corresponding to each frequency band, and perform continuous wavelet transform on each seismic post-stack record. According to the set number of frequency bands and frequency band range, the continuous wavelet transform result is inversely transformed into the time domain to obtain the seismic frequency data required for inversion; Interpolation module: Calculates the fluid mobility coefficient of the reservoir based on the logging curves obtained from well logging and interpolates it into a low-frequency model; Inversion module: Based on relative wave impedance and L2 regularization constraints, the fluid mobility inversion equation is established. The seismic frequency-dividing data, low-frequency model and relative wave impedance are introduced into the established fluid mobility inversion equation, and the mobility inversion results of each seismic post-stack record are solved using gradient descent.

6. The device for frequency-dependent fluid mobility inversion according to claim 5, characterized in that: The specific implementation steps of the earthquake frequency division data acquisition module are as follows: Step 2.1: Set the number of frequency bands to be inverted, and set the frequency range of each frequency band, including the start frequency and end frequency; Step 2.2: Perform continuous wavelet transform on each post-stack earthquake record x(t) corresponding to the frequency band to obtain the continuous wavelet transform result of each post-stack earthquake record x(t). The formula of continuous wavelet transform is: Where W(a, b) represents the continuous wavelet transform result of each post-stack seismic recording signal x(t), a represents the wavelet scale factor, b represents the translation factor, which is used to control the time position, and ψ a，b (t) represents the mother wavelet after translation by b and scaling by a at time t, * represents the complex conjugate, and ψ represents the mother wavelet; Step 2.3: Transform the continuous wavelet transform results of the required frequency band into the time domain to obtain the seismic frequency data required for inversion. The formula is: The corresponding relationship between frequency and scale after continuous wavelet transform is expressed as: W filtered (a,b)=W(a,b)·E(a) (6) Where s(t) is the inverted seismic frequency data, W filtered (a, b) means in a min <a max The new wavelet coefficient result constructed by the corresponding frequency band, C ψ represents the volume constant of the wavelet, f represents the frequency, F c represents the center frequency of the mother wavelet, [f min , f max ] represents the frequency band range, f min <f<f max , a max Indicates the end point of the frequency band f max The corresponding scale, a min Indicates the starting point of the frequency band f min The corresponding scale, E(a) represents the Dirac function.​ 7. The device for frequency-dependent fluid mobility inversion according to claim 6, characterized in that: The specific implementation steps of the interpolation module are: Step 3.1: Obtain the well logging curve based on the well logging curve, based on the P-wave velocity of the reservoir at the i-th sampling point , shear wave velocity Density ρ i and porosity φ i Calculate the fluid mobility coefficient of each sampling point using the following formula: Where C(i) represents the fluid mobility coefficient of the i-th sampling point in the reservoir, R1(i) represents the first-order term of the reflection coefficient of the i-th sampling point in the reservoir, represents the fluid density at the i-th sampling point in the reservoir, Z i represents the impedance of the i-th sampling point in the reservoir, represents the equivalent P-wave velocity of the i-th sampling point in the reservoir, in dry rock M i is the bulk modulus or longitudinal wave modulus of the i-th sampling point in the reservoir, represents the bulk modulus of dry rock at the i-th sampling point in the reservoir, μ i represents the second Lame coefficient of the i-th sampling point in the reservoir, and represents the dimensionless coefficient of the i-th sampling point in the reservoir, represents the propagation velocity of the wave in the fluid at the i-th sampling point in the reservoir, represents the compressibility of the fluid at the i-th sampling point in the reservoir, and denote the elastic moduli of the skeleton and fluid at the i-th sampling point in the reservoir, represents the bulk modulus of the rock matrix at the i-th sampling point in the reservoir, which is obtained using the LRM linear fitting method. i 、D i 、P i , Q i It is an intermediate variable used to simplify the formula during the calculation process; Step 3.2: Interpolate the fluid mobility coefficient of each sampling point in the reservoir into a low-frequency model.

8. The device for frequency-dependent fluid mobility inversion according to claim 7, characterized in that: The specific implementation steps of the inversion module are: Step 4.1: Establish the fluid mobility inversion equation based on relative wave impedance and L2 regularization constraint. The specific steps are as follows: According to Silin, when a fast longitudinal wave is incident vertically, the reflection coefficient at the interface of a viscous fluid-porous medium can be described as follows in the low-frequency asymptotic form: Where R is the longitudinal wave reflection coefficient, ε is a small dimensionless parameter, |ε| = |jλω|, j is an imaginary unit, ω is the frequency, ρ f is the density of viscous fluid, k is the permeability of porous rock, η is the viscosity coefficient of fluid, Z1 and Z2 are the longitudinal wave impedances of the upper and lower media respectively, λ is the intermediate variable, and R1 represents the first-order term of the reflection coefficient; Write equation (20) as a function of time and frequency, and take the real part to get: in, is the mobility coefficient at time t, which means that in the inversion process, the low-frequency model established by calculating the fluid mobility coefficient of all sampling points in the logging reservoir is replaced by K(t). It represents the fluid mobility value of all sampling points of the logging reservoir at time t; Write (21) in matrix form and expand it according to frequency to obtain the following equations: Among them, I n×n is the n-dimensional identity matrix, ω m represents the mth frequency sampling point, is an n-dimensional diagonal matrix; Multiply both ends of equation (22) by the wavelet matrix to obtain: Formula (23) can also be simplified into the following form: G·X0=d (24) Where W is the wavelet matrix, d(t, ω m ) is the earthquake frequency data according to frequency, G represents X0 represents d represents If the seismic wave is incident at a small angle, the relationship between the reflection coefficient and the longitudinal wave impedance is in the form of time and frequency: Where: I P (t, f) is the longitudinal wave impedance in the form of time t and frequency f, I P (t0, f) is the initial impedance value, Combining equations (22) and (25) and expanding them according to frequency, we get the model constraint equation: in, It is expanded into frequency-dependent impedance, D is the summation matrix, and Equation (26) is simplified into the following form: C′·X0=ε (27) Where C′ represents X0 represents ε represents Combining equations (27) and (24), we construct the inversion objective function: Introducing the L2 norm into the inversion objective function yields a new objective function: Where α is the impedance model Parameters, β is the regularization coefficient; Step 4.2: Gradient descent method is used to solve the fluid mobility, that is, to solve equation (29) based on the gradient descent method. The specific steps are as follows: Based on the new objective function (29), the gradient formula is: Where G T represents the transpose of the G matrix, C′ T It is represented as the transpose of the C′ matrix, and E is represented as the identity matrix; The gradient of the first round is obtained based on the gradient formula: If satisfied Then the mobility inversion results of each earthquake post-stack record are obtained Otherwise, the update using gradient descent is Then make a judgment; Where η is the learning rate; this formula indicates that when the second norm of the gradient is less than the threshold ε, the gradient descent algorithm stops, indicating that it is close to the optimal solution.

Images