WO2021147529A1 - 基于更新概率比率恒定理论的多点地质统计叠前反演方法 - Google Patents

基于更新概率比率恒定理论的多点地质统计叠前反演方法 Download PDF

Info

Publication number
WO2021147529A1
WO2021147529A1 PCT/CN2020/134088 CN2020134088W WO2021147529A1 WO 2021147529 A1 WO2021147529 A1 WO 2021147529A1 CN 2020134088 W CN2020134088 W CN 2020134088W WO 2021147529 A1 WO2021147529 A1 WO 2021147529A1
Authority
WO
WIPO (PCT)
Prior art keywords
data
lithofacies
grid
estimated
point
Prior art date
Application number
PCT/CN2020/134088
Other languages
English (en)
French (fr)
Inventor
尹艳树
胡迅
冯文杰
张昌民
Original Assignee
长江大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 长江大学 filed Critical 长江大学
Priority to US17/315,503 priority Critical patent/US11131784B2/en
Publication of WO2021147529A1 publication Critical patent/WO2021147529A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B47/00—Survey of boreholes or wells
    • E21B47/002—Survey of boreholes or wells by visual inspection
    • E21B47/0025—Survey of boreholes or wells by visual inspection generating an image of the borehole wall using down-hole measurements, e.g. acoustic or electric
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • E—FIXED CONSTRUCTIONS
    • E21—EARTH OR ROCK DRILLING; MINING
    • E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
    • E21B2200/00—Special features related to earth drilling for obtaining oil, gas or water
    • E21B2200/20—Computer models or simulations, e.g. for reservoirs under production, drill bits
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6169—Data from specific type of measurement using well-logging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters

Definitions

  • the embodiments of the present disclosure relate to the technical field of oil and gas exploration and development, and more specifically, a multi-point geostatistical pre-stack inversion method based on the theory of constant update probability ratio.
  • Seismic random inversion is the fusion of logging data under the constraints of seismic data, which can broaden the frequency spectrum of the inversion results, improve the resolution, and search for thinner and smaller reservoirs.
  • Hass and Dubrule (1994) first proposed the prototype of random inversion, combining the sequential Gaussian simulation method with seismic inversion.
  • Yin Xingyao et al. (2005) applied Bayesian theory to inversion and got good results.
  • Sun Yuecheng applied the stochastic inversion method to the pre-stack inversion.
  • Seismic random inversion methods have gradually matured in the past two decades of research, but there are still many shortcomings, such as low computational efficiency and poor lateral continuity. Therefore, more in-depth research on random inversion methods is needed. , So as to better predict the reservoir and improve the accuracy of exploration and development.
  • the purpose of the embodiments of the present disclosure is to overcome the shortcomings of the above-mentioned background technology and propose a multi-point geostatistical pre-stack inversion method based on the theory of constant update probability ratio.
  • the seismic data is pre-stack seismic data
  • the well data includes drilling core and logging data
  • the lithofacies corresponding to different well depths and the elastic parameters of different lithofacies can be determined
  • the elastic parameters include density, longitudinal wave velocity and shear wave velocity; then establish the comparison of different lithofacies according to the lithofacies data, and establish the cumulative distribution function of different lithofacies elastic parameters according to the elastic parameter data; establish training images that meet the characteristics of the reservoir in the work area;
  • Core data and logging data are allocated to the nearest grid node as hard data
  • the information contained in the actual drilling within the scope of the work area is called hard data or conditional data.
  • the hard data includes lithofacies, density, longitudinal wave data, and shear wave velocity; the grid with hard data allocated to lithofacies is called Known grids or known points; grids without lithofacies data are called grids to be estimated or points to be estimated;
  • the allocated well data is drilling data and logging data
  • the grid is a grid allocated with lithofacies data
  • the known grid is the grid through which the well is drilled, and the grid to be estimated is the grid that has not been penetrated by the well;
  • the known grid in the work area contains lithofacies, density, P-wave velocity and shear wave velocity, and the grid to be estimated does not have any data;
  • initial attribute values are assigned to the points to be estimated in the work area
  • the shape of the data template can be determined according to the heterogeneity of the sedimentary facies.
  • a two-dimensional ellipse or a three-dimensional ellipsoid can be used.
  • Two-dimensional rectangles or three-dimensional cuboids can be used when the performance is weak; there is no hard and fast rule on the size of the data template, and the data template size of 5 ⁇ 7 or 5 ⁇ 5 ⁇ 7 can generally be used;
  • the inversion is obtained through the following 3 sub-steps:
  • Sub-step 5.1 Visit the simulation node sequentially
  • Sub-step 5.2 Suggested data model acquisition
  • the initial inversion when processing the grid to be estimated, you must first determine the number of conditional points within the data template centered on the point to be estimated, the relative position of the point to be estimated, and the lithofacies type to form the point to be estimated Data events presented in vector form as the center; use data events to randomly scan the training images established in step 1, and obtain the first completely matching data pattern as the suggested data pattern; the suggested data pattern represents a deposition pattern, which has data The lithofacies data structure of the sample size, the lithofacies data structure of the condition points in the data template can be fully reflected in the recommended data model;
  • a data event when processing the grid to be estimated, a data event must be formed according to the point to be estimated and the conditional point of the data template range of the point to be estimated, and the training image established in step 1 is randomly scanned using the data event to obtain several complete data.
  • the matched data pattern is used as the candidate data pattern; the candidate data patterns are sorted according to the lithofacies type of the center point to make a cumulative probability distribution;
  • the candidate data pattern can include at most the first 50 data patterns scanned from the first-to-last training image.
  • all data patterns of the training image are scanned by the data event as the candidate data pattern, according to
  • the different lithofacies in the center of the candidate data model calculates the ratio of different lithofacies, which is recorded as a phase 1 : a phase 2 : ...: a phase m , and then calculate the probability of the elastic parameters of different phases according to the following formula according to the last elastic parameter:
  • b represents the probability of the lithofacies elastic parameter
  • n represents the number of grids in the data template range of the point to be estimated
  • ⁇ i represents the density value of the i-th grid in the range of the data template around the point to be estimated.
  • phase 1 ⁇ b phase 1 (a phase 2 ⁇ b phase 2 ): ...: (a phase m ⁇ b phase m ), denoted as P phase 1 : P phase 2 : ...: P facies m ; normalize the joint probability to obtain P'facies 1 : P'facies 2 : ...: P'facies m , and make the joint probabilities of different lithofacies into a cumulative distribution function;
  • Sub-step 5.3 Suggested acquisition of elastic parameters
  • the rock facies at each position in the recommended data model is sampled for elastic parameters;
  • the elastic parameters include density values, longitudinal wave velocity values, and shear wave velocity values;
  • Sub-step 5.4 Calculate the composite record of the recommended elasticity parameter in the recommended data model
  • EI( ⁇ ) represents the reflection coefficient when the incident angle is ⁇
  • ⁇ represents the preset incident angle
  • V p represents the longitudinal wave velocity
  • V s represents the shear wave velocity
  • ⁇ represents the density
  • g(m) represents the initial forward model
  • m represents the elastic parameter
  • w represents the seismic wavelet
  • Sub-step 5.5 Select the elastic parameter with the highest matching rate with the original pre-stack data
  • the last elasticity parameter and the recommended elasticity parameter are calculated in steps 5.3 and 5.4 respectively to calculate the composite record, the former is recorded as g (m last time ), the latter is recorded as g (m suggestion ); here the last elastic parameter can be the initial elasticity
  • the parameter may be the suggested elasticity parameter during the previous iteration;
  • m represents the selected elastic parameter
  • m last represents the last elastic parameter
  • m suggestion represents the recommended elastic parameter in the recommended data model
  • step 5 When the number of iterative inversions is less than 7, repeat step 5 and perform the next iterative inversion; when the iterative inversion is completed for the seventh time, the entire inversion process ends, and a multi-point geostatistical seismic inversion model is obtained. Output simulation results.
  • step 2 when there are multiple well data in a grid, the data point closest to the grid center is selected as the hard data.
  • step 5.1 the order of accessing the simulation nodes is usually from the area with rich well data to the area with less well data, and finally to the area without wells.
  • the rock facies condition data includes well condition data and sedimentary facies data of simulated nodes.
  • step 5.2 when it is not the first iteration, it is recommended that the data mode be jointly determined based on the joint probability obtained from the prior geological information of multi-point geostatistics and the cumulative distribution function of the elastic parameters of the logging data.
  • the initial attribute value assigned to the point to be estimated in the work area can be a definite value, or it can be randomly obtained by Monte Carlo sampling.
  • step 1 the cumulative distribution function is used as the object of the elastic parameter sampling under the constraints of the later lithofacies.
  • step 3 in step 3;
  • the initial attribute value assigned to the point to be estimated in the work area is obtained through the following three sub-steps:
  • Sub-step 3.1 Sequential access to simulation nodes
  • Visit the simulation node to traverse all grids according to a certain order, such as from bottom to top or from left to right or from front to back; when there is lithofacies data in the grid, there is no need to assign initial attribute values; in the grid When there is no lithofacies data, proceed to the next step;
  • Sub-step 3.2 Randomly allocate lithofacies according to Monte Carlo sampling
  • the Monte Carlo sampling method is used to randomly sample to allocate lithofacies to the estimated grid.
  • the lithofacies are only temporarily stored in the grid, and the elasticity is allocated according to the lithofacies in the next step. After the parameters can be removed;
  • Sub-step 3.3 Randomly assign the corresponding elastic parameters of rocks to different lithofacies according to Monte Carlo sampling;
  • Monte Carlo sampling is performed to obtain the density value, longitudinal wave velocity and value of the grid temporary lithofacies. Shear wave velocity value;
  • the known grids of the work area are still the grids that the wells have passed through, and the grids to be estimated are still the grids that have not been drilled through; and the lithofacies data is only available on the known grids.
  • the elastic parameter data is available in all grids in the work area.
  • the embodiments of the present disclosure have the following advantages: 1.
  • the embodiments of the present disclosure use a multi-point geostatistical method to obtain prior information, and then screen through a minimum objective function, thereby reducing the complexity of seismic inversion.
  • a single fully matched suggested data pattern simulation obtained by scanning a training image through a data event allows the calculation efficiency to be greatly improved on the premise that the accuracy of the inversion result is ensured.
  • Fig. 1 is a step diagram of Embodiment 1 of an embodiment of the disclosure.
  • Figure 2 is the well condition data of Example 1 of the embodiments of the disclosure.
  • Fig. 3 is a density cumulative distribution diagram of Example 1 of an embodiment of the disclosure.
  • FIG. 4 is a cumulative distribution diagram of longitudinal wave velocity in Embodiment 1 of the embodiment of the disclosure.
  • FIG. 5 is a cumulative distribution diagram of transverse wave velocity in Embodiment 1 of the embodiment of the disclosure.
  • FIG. 6 is a training image of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 7 is a data template of Embodiment 1 of the embodiments of the disclosure.
  • FIG. 8 is a data event in Embodiment 1 of the embodiment of the disclosure.
  • FIG. 9 is a suggested data pattern of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 10 is a partial diagram of the pseudo-random number in Embodiment 1 of the embodiment of the disclosure.
  • FIG. 11 is the first iterative inversion phase diagram of Embodiment 1 of the embodiments of the disclosure.
  • FIG. 12 is the first iterative inversion density map of Embodiment 1 of the embodiments of the disclosure.
  • FIG. 13 is the first iterative inversion of longitudinal wave velocity in Embodiment 1 of the embodiments of the disclosure.
  • FIG. 14 is the first iterative inversion of the shear wave velocity diagram in Embodiment 1 of the embodiments of the disclosure.
  • FIG. 15 is a first iterative inversion of synthetic seismic records in Embodiment 1 of the disclosed embodiments.
  • FIG. 16 is a second iteration inversion phase diagram of Embodiment 1 of the embodiments of the disclosure.
  • FIG. 17 is a second iterative inversion density map of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 18 is a diagram of the longitudinal wave velocity inversion in the second iteration of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 19 is a diagram of the transverse wave velocity inversion for the second iteration in Embodiment 1 of the embodiment of the disclosure.
  • FIG. 20 is a second iterative inversion of the synthetic seismic record in Embodiment 1 of the embodiment of the disclosure.
  • Fig. 21 is a seventh iterative inversion phase diagram of Embodiment 1 of the embodiments of the disclosure.
  • FIG. 22 is a seventh iterative inversion density map of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 23 is a graph of longitudinal wave velocity inversion in the seventh iteration of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 24 is a graph of the transverse wave velocity inversion in the seventh iteration of Embodiment 1 of the embodiment of the disclosure.
  • FIG. 25 is a seventh iterative inversion of synthetic seismic record in Embodiment 1 of the embodiment of the disclosure.
  • the multi-point geostatistical pre-stack inversion method based on the theory of constant update probability ratio described in the embodiment of the present disclosure includes obtaining training images of the reservoir structure in the study area through geological analysis of field outcrops, geological knowledge base, well data, and seismic data; The data is hard-constrained, and the multi-point modeling method is used to obtain the prior information of the lithofacies distribution; through the probability density distribution relationship between the attributes and the lithofacies, Monte Carlo sampling obtains the seismic attributes of the predicted points; when it is not the first iteration, the previous iteration
  • the stratigraphic model parameters are used as constraints, and the multi-point modeling method uses joint probability to obtain the probability density distribution of lithofacies;
  • the seismic data is pre-stack seismic data
  • the well data includes drilling core and logging data
  • the lithofacies corresponding to different well depths and the elastic parameters of different lithofacies can be determined
  • the elastic parameters include density, longitudinal wave velocity and shear wave velocity; then establish the comparison of different lithofacies according to the lithofacies data, and establish the cumulative distribution function of different lithofacies elastic parameters according to the elastic parameter data; establish training images that meet the characteristics of the reservoir in the work area;
  • Core data and logging data are allocated to the nearest grid node as hard data
  • the information contained in the actual drilling within the scope of the work area is called hard data or conditional data.
  • the hard data includes lithofacies, density, longitudinal wave data, and shear wave velocity; the grid with hard data allocated to lithofacies is called Known grids or known points; grids without lithofacies data are called grids to be estimated or points to be estimated;
  • the allocated well data is drilling data and logging data
  • the grid is a grid allocated with lithofacies data
  • the known grid is the grid through which the well is drilled, and the grid to be estimated is the grid that has not been penetrated by the well;
  • the known grid in the work area contains lithofacies, density, P-wave velocity and shear wave velocity, and the grid to be estimated does not have any data;
  • initial attribute values are assigned to the points to be estimated in the work area
  • the shape of the data template can be determined according to the heterogeneity of the sedimentary facies.
  • a two-dimensional ellipse or a three-dimensional ellipsoid can be used.
  • Two-dimensional rectangles or three-dimensional cuboids can be used when the performance is weak; there is no hard and fast rule on the size of the data template, and the data template size of 5 ⁇ 7 or 5 ⁇ 5 ⁇ 7 can generally be used;
  • the inversion is obtained through the following 3 sub-steps:
  • Sub-step 5.1 Visit the simulation node sequentially
  • Sub-step 5.2 Suggested data model acquisition
  • the initial inversion when processing the grid to be estimated, you must first determine the number of conditional points within the data template centered on the point to be estimated, the relative position of the point to be estimated, and the lithofacies type to form the point to be estimated Data events presented in vector form as the center; use the data events to randomly scan the training images established in step 1, and obtain the first completely matched data pattern as the recommended data pattern; the recommended data pattern represents a deposition pattern, which is a data pattern
  • the lithofacies data structure of the sample size, the lithofacies data structure of the condition points in the data template can be fully reflected in the recommended data model;
  • a data event when processing the grid to be estimated, a data event must be formed according to the point to be estimated and the conditional point of the data template range of the point to be estimated, and the training image established in step 1 is randomly scanned using the data event to obtain several complete data.
  • the matched data pattern is used as the candidate data pattern; the candidate data patterns are sorted according to the lithofacies type of the center point to make a cumulative probability distribution;
  • the candidate data pattern can include at most the first 50 data patterns scanned from the first-to-last training image.
  • all data patterns of the training image are scanned by the data event as the candidate data pattern, according to
  • the different lithofacies in the center of the candidate data model calculates the ratio of different lithofacies, which is recorded as a phase 1 : a phase 2 : ...: a phase m , and then calculate the probability of the elastic parameters of different phases according to the following formula according to the last elastic parameter:
  • b represents the probability of the lithofacies elastic parameter
  • n represents the number of grids in the data template range of the point to be estimated
  • ⁇ i represents the density value of the i-th grid in the range of the data template around the point to be estimated.
  • phase 1 ⁇ b phase 1 (a phase 2 ⁇ b phase 2 ): ...: (a phase m ⁇ b phase m ), denoted as P phase 1 : P phase 2 : ...: P facies m ; normalize the joint probability to obtain P'facies 1 : P'facies 2 : ...: P'facies m , and make the joint probabilities of different lithofacies into a cumulative distribution function;
  • Sub-step 5.3 Suggested acquisition of elastic parameters
  • the rock facies at each position in the recommended data model is sampled for elastic parameters;
  • the elastic parameters include density values, longitudinal wave velocity values, and shear wave velocity values;
  • Sub-step 5.4 Calculate the composite record of the recommended elasticity parameter in the recommended data model
  • EI( ⁇ ) represents the reflection coefficient when the incident angle is ⁇
  • ⁇ represents the preset incident angle
  • V p represents the longitudinal wave velocity
  • V s represents the shear wave velocity
  • ⁇ represents the density
  • g(m) represents the initial forward model
  • m represents the elastic parameter
  • w represents the seismic wavelet
  • Sub-step 5.5 Select the elastic parameter with the highest matching rate with the original pre-stack data
  • the last elasticity parameter and the recommended elasticity parameter are calculated in steps 5.3 and 5.4 respectively to calculate the composite record, the former is recorded as g (m last time ), the latter is recorded as g (m suggestion ); here the last elastic parameter can be the initial elasticity
  • the parameter may be the suggested elasticity parameter during the previous iteration;
  • m represents the selected elastic parameter
  • m last represents the last elastic parameter
  • m suggestion represents the recommended elastic parameter in the recommended data model
  • step 5 When the number of iterative inversions is less than 7, repeat step 5 and perform the next iterative inversion; when the iterative inversion is completed for the seventh time, the entire inversion process ends, and a multi-point geostatistical seismic inversion model is obtained. Output simulation results.
  • step 2 when there are multiple well data in a grid, the data point closest to the grid center is selected as the hard data.
  • the order of accessing the simulation nodes is usually from the area with rich well data to the area with less well data, and finally to the area without wells.
  • the rock facies condition data includes well condition data and sedimentary facies data of simulated nodes.
  • sub-step 5.2 when it is not the first iteration, it is recommended that the data model be jointly determined based on the joint probability obtained from the prior geological information of multi-point geostatistics and the cumulative distribution function of the elastic parameters of the logging data.
  • the initial attribute value assigned to the point to be estimated in the work area can be a definite value, or it can be randomly obtained by Monte Carlo sampling.
  • step 1 the cumulative distribution function is used as the object of elastic parameter sampling under the constraints of later lithofacies.
  • step 3 the initial attribute value of the point to be estimated in the work area is obtained through the following three sub-steps:
  • Sub-step 3.1 Sequential access to simulation nodes
  • Visit the simulation node to traverse all grids according to a certain order, such as from bottom to top or from left to right or from front to back; when there is lithofacies data in the grid, there is no need to assign initial attribute values; in the grid When there is no lithofacies data, proceed to the next step;
  • Sub-step 3.2 Randomly allocate lithofacies according to Monte Carlo sampling
  • the Monte Carlo sampling method is used to randomly sample to allocate lithofacies to the estimated grid.
  • the lithofacies are only temporarily stored in the grid, and the elasticity is allocated according to the lithofacies in the next step. After the parameters can be removed;
  • Sub-step 3.3 Randomly assign the corresponding elastic parameters of rocks to different lithofacies according to Monte Carlo sampling;
  • Monte Carlo sampling is performed to obtain the density value, longitudinal wave velocity and value of the grid temporary lithofacies. Shear wave velocity value;
  • the known grids of the work area are still the grids that the wells have passed through, and the grids to be estimated are still the grids that have not been drilled through; and the lithofacies data is only available on the known grids.
  • the elastic parameter data is available in all grids in the work area.
  • Embodiment 1 Refer to Figure 1 for the implementation step diagram of the seismic reservoir inversion method described in this embodiment, specifically:
  • the training image shown in Figure 2 is established according to the features of the river facies in the work area. Establish a cumulative distribution map of the elastic parameters of sand and mudstone rocks in the work area.
  • the density range of mudstone is 2.17-2.33g/cm3
  • the longitudinal wave velocity range is 3930-4272m/s
  • the shear wave velocity range is 2226-2496m/s
  • the density range of sandstone is 2.27- 2.42g/cm3
  • the longitudinal velocity range is 4041-4528m/s
  • the lateral velocity range is 2300-2561m/s.
  • the physical parameters of sand and mudstone overlap and each obey Gaussian distribution.
  • the total length of the work area is 1000m, and the total thickness is 45m.
  • Divide the model grid into 200 ⁇ 1 ⁇ 45, each of which has a width of 5 meters and a thickness of 1m.
  • the facies data and rock elastic parameters of the actual logging are allocated to the grid of the work area.
  • the initial attribute value is assigned to the point to be estimated in the work area, and the initial value is extracted from the statistical probability density distribution using the Monte Carlo method.
  • the phase diagram, density diagram, P-wave velocity diagram, shear-wave velocity diagram and synthetic seismic record of the second iteration inversion are shown in Figures 14, 15, 16, 17 and 18.
  • the difference from the first iterative inversion is that after extracting the first 50 data patterns as data patterns, they are sorted according to the relative data patterns of the center points to form a cumulative distribution map. For example, there are 15 data patterns where the center grid is mudstone facies. , There are 35 data models for sandstone facies.
  • the grid density value is 2.32 after the last iteration, the longitudinal wave velocity is 4320, and the shear wave velocity is 2461, then according to the probability density function of sand and mudstone, the ratio of mudstone facies to sandstone facies is 1:1, so many points and The joint probability ratio of earthquakes is 3:7.
  • Monte Carlo sampling is performed to determine whether the central grid is a data model of mudstone or sandstone facies. For example, when the sampling result is 0.4, the cumulative distribution map of the selected data model corresponds to the 20th one The data model centered on sandstone facies is used as the recommended data model.
  • the elastic parameters of the lithofacies at different positions of the proposed data model are sampled, the reflection coefficient is calculated according to the Connolly formula, and the reflection coefficient and seismic wavelet convolution are calculated to synthesize seismic records. Compare the synthetic seismic record with the original seismic record and select the best mode and its corresponding attribute retention. Simulate all grids in turn.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Remote Sensing (AREA)
  • Geology (AREA)
  • Environmental & Geological Engineering (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • General Physics & Mathematics (AREA)
  • Mining & Mineral Resources (AREA)
  • Fluid Mechanics (AREA)
  • Geochemistry & Mineralogy (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种基于更新概率比率恒定理论的多点地质统计叠前反演方法，包括如下步骤：整理资料、工区网格化和分配井数据、赋予模拟工区初始属性值、选择适当大小的数据样板、反演和判断迭代终止；通过利用多点地质统计学方法获得先验信息，然后通过最小目标函数进行筛选，降低了地震反演的复杂性，有效提高了储层反演的精度。

Description

基于更新概率比率恒定理论的多点地质统计叠前反演方法 技术领域
本公开实施例涉及到油气勘探开发技术领域，更加具体地是基于更新概率比率恒定理论的多点地质统计叠前反演方法。
背景技术
如今，油气勘探开发的难度越来越大，勘探程度越来越高，寻找新的油气藏越来越难，这就要求储层参数得到的要准确，使勘探的不确定性降低一些。
地震随机反演是在地震数据的约束下，融合测井数据，可以拓宽反演结果的频谱，提高分辨率，一遍于寻找更薄更小的储层。Hass和Dubrule(1994)最早提出了随机反演的雏形，将序贯高斯模拟方法结合地震反演。印兴耀等人(2005)将贝叶斯理论应用于反演中，得到了不错的结果。2010年孙月成将随机反演方法应用到叠前反演中。2011年王家华研究了地震数据的多点统计学方法，他将井数据作为硬数据，地震数据作为约束条件，约束随机模拟过程，降低了井间的不确定性。2011年尹艳树总结了多点地质统计学方法的应用。
地震随机反演方法在近二十年的研究中逐渐走向成熟，但是还有很多缺陷，比如计算效率低、横向连续性差等问题，因此，还需要不断地对随机反演方法进行更加深入的研究，从而更好的进行储层的预测和提高勘探开发的精度。
因此，迫切需要一种进行储层的预测和提高勘探开发的精度的方法。
发明内容
本公开实施例的目的在于克服上述背景技术的不足之处，而提出基于更新概率比率恒定理论的多点地质统计叠前反演方法。
本公开实施例的目的是通过如下技术方案来实施的：基于更新概率比率恒定理论的多点地质统计叠前反演方法，它包括如下具体步骤；
①、整理资料
检查原始地震资料和井资料否齐全，其中地震资料为叠前地震数据，井资料包括钻井岩心和测井数据；根据井资料解释成果，可以确定不同井深对应的岩相以及不同岩相的弹性参数，弹性参数包括密度、纵波速度和横波速度；接着根据岩相数据建立不同岩相的相比例，根据弹性参数数据建立不同岩相弹性参数的累积分布函数；建立符合工区储层特征的训练图像；
②、工区网格化和分配井数据
根据实际工区范围选择合适的网格尺寸大小，对工区进行网格划分，建立网格模型；根据各井的平面位置和井深将岩心数据和测井数据网格化；岩心数据即为岩相，测井数据为解释的密度、纵波速度和横波速度；
岩心数据和测井数据作为硬数据分配到划分的最邻近的网格节点上；
将工区范围内实际钻井上包含的信息称为硬数据或者称为条件数据，所述的硬数据包括了岩相、密度、纵波数据和横波速度；将分配了岩相的硬数据的网格称为已知网格或者称已知点；没有岩相数据的网格称为待估网格或又叫待估点；
在工区网格化和分配井数据后，所述的分配井数据为钻井数据和测井数据；
已知网格为分配有岩相数据的网格，
已知网格就是钻井穿过的网格，待估网格就是没有钻井穿过的网格；工区已知网格包含岩相、密度、纵波速度以及横波速度，待估网格没有任何数据；
③、赋予模拟工区初始属性值
根据统计的工区岩石弹性参数累积分布函数，对工区待估点赋予初始属性值(指弹性参数)；
④、选择适当大小的数据样板
根据沉积相形态特征确定数据样板形状和大小；数据样板的形状可以根据沉积相的非均质性来定，非均质性强时可以采用二维的椭圆形或者三维的椭球状，非均质性较弱时可以采用二维长方形或者三维长方体；数据样板的大小没有硬性规定，一般可以采用5×7或者是5×5×7的数据样板大小；
⑤、反演
反演通过下面3个分步骤得到：
分步骤5.1：顺序访问模拟节点
先对工区所有网格随机产生伪随机数，该随机数的数值在0到1之间，然后以待估网格的数据样板大小为周围搜索条件点个数；每个待估点的数值等于伪随机数与条件点个数之和，对所有的待估点其进行从大到小排序，值较大的网格优先模拟，从而得到整个工区的模拟路径；
分步骤5.2：建议数据模式获取
如果进行初次反演，处理待估网格时，先要确定以待估点为中心的数据样板范围内的条件点的个数、相对待估点的位置以及岩相类型，形成以待估点为中心的以矢量形式呈现的数据事件；利用数据事件随机扫描步骤①中建立的训练图像，从中获得首个完全匹配的数据模式作为建议数据模式；建议数据模式表示一种沉积模式，是具有数据样板大小的岩相数据结构，数据样板中条件点的岩相数据结构在建议数据模式中可以完全体现；
如果进行迭代反演，处理待估网格时，先要根据待估点和待估点数据样板范围的条件点形成数据事件，利用数据事件随机扫描步骤①中建立的训练图像，从中获得若干完全匹配的数据模式作为候选数据模式；将候选数据模式按照中心点岩相类型排序，做成累积概率分布；
候选数据模式最多可包含从先到后训练图像中扫描到的前50个数据模式，当不到50个候选数据模式时，则用数据事件扫描完训练图像所有的数据模式作为候选数据模式，根据候选数据模式中心岩相的不同计算不同岩相的比例，记做a 相1：a 相2：…：a 相m，再根据上一次弹性参数根据下述公式计算不同相的弹性参数的概率：
Figure PCTCN2020134088-appb-000001
公式中，b表示岩相弹性参数的概率，n表示待估点数据样板范围网格个数，ρ i表示待估点周围数据样板范围内第i个网格的密度值，
Figure PCTCN2020134088-appb-000002
表示待估点周围数据样板范围内第i个网格的纵波速度，
Figure PCTCN2020134088-appb-000003
表示待估点周围数据样板范围内第i个网格的横波速度，μ ρ表示该岩相密度均值，μ Vp表示该岩相纵波速度均值，μ Vs表示该岩相横波速度均值，
Figure PCTCN2020134088-appb-000004
表示该岩相密度方差，
Figure PCTCN2020134088-appb-000005
表示该岩相纵波速度方差，
Figure PCTCN2020134088-appb-000006
表示该岩相横波速度方差；
不同岩相的弹性参数的比例记做b 相1：b 相2：…：b 相m；
最后不同相的联合概率用(a 相1·b 相1)：(a 相2·b 相2):…:(a 相m·b 相m)计算，记做P 相1：P 相2：…：P 相m；将联合概率归一化处理，得到P’ 相1：P’ 相2：…：P’ 相m，将不同岩相的联合概率做成累积分布函数；
接着将候选数据模式按照中心网格的岩相进行排序，再随机抽样获取0到1之间的随机数，根据随机数在联合概率形成累积分布函数的位置和候选模式的累积概率分布，来确定建议数据模式；
分步骤5.3：建议弹性参数的获取
根据步骤①统计的不同岩石相弹性参数的累积分布函数，对建议数据模式中的各位置的岩相进行弹性参数的抽样；所述的弹性参数包括密度值和纵波速度值和横波速度值；
分步骤5.4：计算建议数据模式中建议弹性参数的合成记录
将建议弹性参数与预设入射角度基于下述Connolly公式计算反射系数；
Figure PCTCN2020134088-appb-000007
公式中，
Figure PCTCN2020134088-appb-000008
EI(θ)表示入射角为θ时的反射系数，θ代表预设入射角度，V p代表纵波速度，V s代表横波速度，ρ代表密度；
基于反射系数和地震地波，采用下述褶积公式计算初始正演模型：
g(m)＝w*EI(θ)
公式中，g(m)表示初始正演模型，m表示弹性参数，w表示地震子波；
分步骤5.5：选择和原始叠前数据匹配率最高的弹性参数
将上一次弹性参数和建议弹性参数通过分步骤5.3和5.4分别计算合成记录，前者记为g(m 上次)，后者记为g(m 建议)；这里上一次的弹性参数可以是初始弹性参数或可以是前一次迭代过程中的建议弹性参数；
采用如下公式将若干条合成地震道与实际地震记录比较，选择目标函数最小的地层模型参数和对应的建议数据模式作为最优反演结果：
Figure PCTCN2020134088-appb-000009
公式中，m代表选择的弹性参数，m 上次表示上一次的弹性参数，m 建议表示建议数据模型中建议弹性参数；
⑥、判断迭代终止
重复分步骤5.1-分步骤5.5，直到遍历所有网格，即所有待估点处模拟完成，就是实现了一次反演实现；
当迭代反演次数小于7次时，则重复步骤⑤，进行下一次迭代反演；当迭代反演第7次完成，则整个反演过程结束，得到多点地质统计学地震反演的模型，输出模拟结果。
在上述技术方案中：在步骤②中；当一个网格中存在多个井数据时，选择离网格中心最近的数据点作为硬数据。
在上述技术方案中：在分步骤5.1中；访问模拟节点的顺序通常是由井资料丰富区域逐步过渡到井资料较少区域，最后到无井区域。
在上述技术方案中：在分步骤5.2中；岩石相条件数据包括井条件数据和已模拟节点的沉积相数据。
在上述技术方案中：在分步骤5.2中；当不是初次迭代时，建议数据模式根据多点地质统计的先验地质信息和测井资料的弹性参数累积分布函数得到的联合概率共同决定。
在上述技术方案中：在步骤③中；工区待估点赋予初始属性值可以是确定值，也可以是蒙特卡洛抽样随机获取。
在上述技术方案中：在步骤①中；累积分布函数作为后期岩相约束下弹性参数抽样的对象。
在上述技术方案中：在步骤③中；
对工区待估点赋予初始属性值通过下面三个分步骤得到：
分步骤3.1：顺序访问模拟节点
访问模拟节点根据从一定顺序遍历所有网格，比如从下到上或从左到右或从前到后的原则；在网格中有岩相数据时，就不用赋予初始属性值；在网格中没有岩相数据时，则进行下一步；
分步骤3.2：根据蒙特卡洛抽样随机分配岩相；
根据步骤①中统计的不同岩相比例，采用蒙特卡洛抽样方法随机抽样，来对待估网格分配岩相，此时岩相只是临时保存在网格中，在下一分步骤根据岩相分配弹性参数后即可去掉；
分步骤3.3：根据蒙特卡洛抽样随机对不同岩相随机分配岩相对应的弹性参数；
根据分步骤3.2中网格的临时岩相数据，根据步骤①的不同岩相的不同弹性参数数据的累积概率分布，进行蒙特卡洛抽样获得该网格临时岩相的密度值、纵波速度值和横波速度值；
在赋予模拟工区初始属性值后，工区已知网格仍然是钻井穿过的网格，待估网格仍然是没有钻井穿过的网格；而岩相数据仅在已知网格上有，而弹性参数数据在工区所有网格中都有。
本公开实施例具有如下优点：1、本公开实施例利用多点地质统计学方法获得先验信息，然后通过最小目标函数进行筛选，降低了地震反演的复杂性。
2、本公开实施例通过数据事件扫描训练图像获得的单个完全匹配的建议数据模式模拟，使得在确保了反演结果准确性的前提下，大大提高运算效率。
3、本公开实施例由于每次更新均是同时完成了沉积相与弹性参数的更新，因此较好的保证了储层反演一致性；且由于每一次内循环中选取的目标函数最小，即合成地震数据更接近真实地震记录，有效提高了储层反演的精度。
附图说明
图1为本公开实施例实施例1的步骤图。
图2为本公开实施例实施例1的井条件数据。
图3为本公开实施例实施例1的密度累积分布图。
图4为本公开实施例实施例1的纵波速度累积分布图。
图5为本公开实施例实施例1的横波速度累积分布图。
图6为本公开实施例实施例1的训练图像。
图7为本公开实施例实施例1的数据样板。
图8为本公开实施例实施例1的一个数据事件。
图9为本公开实施例实施例1的建议数据模式。
图10为本公开实施例实施例1的伪随机数的局部图。
图11为本公开实施例实施例1第一次迭代反演相图。
图12为本公开实施例实施例1第一次迭代反演密度图。
图13为本公开实施例实施例1第一次迭代反演纵波速度图。
图14为本公开实施例实施例1第一次迭代反演横波速度图。
图15为本公开实施例实施例1第一次迭代反演合成地震记录图。
图16为本公开实施例实施例1第二次迭代反演相图。
图17为本公开实施例实施例1第二次迭代反演密度图。
图18为本公开实施例实施例1第二次迭代反演纵波速度图。
图19为本公开实施例实施例1第二次迭代反演横波速度图。
图20为本公开实施例实施例1第二次迭代反演合成地震记录图。
图21为本公开实施例实施例1第七次迭代反演相图。
图22为本公开实施例实施例1第七次迭代反演密度图。
图23为本公开实施例实施例1第七次迭代反演纵波速度图。
图24为本公开实施例实施例1第七次迭代反演横波速度图。
图25为本公开实施例实施例1第七次迭代反演合成地震记录图。
具体实施方式
以下结合附图详细说明本公开实施例的实施情况，但它们不构成对本公开实施例的限定，仅作举例而已。同时通过和具体实施对本公开实施例作进一步的详细描述。
本公开实施例所述基于更新概率比率恒定理论的多点地质统计叠前反演方法包括通过野外露头、地质知识库、井资料以及地震资料地质分析，获得研究区储层结构训练图像；以井数据为硬约束条件，利用多点建模方法获取岩相分布先验信息；通过属性与岩相的概率密度分布关系，蒙特卡洛抽样获得预测点的地震属性；非初次迭代时，以上一次迭代的地层模型参数作为约束条件，和多点建模方法用联合概率来获取岩相的概率密度分布；
根据Connolly公式计算各道反射系数，并与实际地震子波褶积合成地震道；利用目标函数确定合成记录和实际地震记录的匹配率；最终选取匹配率高的地层模型参数，来实现基于更新概率比率恒定理论的多点地质统计叠 前反演方法。
参照图1-25所示：基于更新概率比率恒定理论的多点地质统计叠前反演方法，它包括如下具体步骤；
①、整理资料
检查原始地震资料和井资料否齐全，其中地震资料为叠前地震数据，井资料包括钻井岩心和测井数据；根据井资料解释成果，可以确定不同井深对应的岩相以及不同岩相的弹性参数，弹性参数包括密度、纵波速度和横波速度；接着根据岩相数据建立不同岩相的相比例，根据弹性参数数据建立不同岩相弹性参数的累积分布函数；建立符合工区储层特征的训练图像；
②、工区网格化和分配井数据
根据实际工区范围选择合适的网格尺寸大小，对工区进行网格划分，建立网格模型；根据各井的平面位置和井深将岩心数据和测井数据网格化；岩心数据即为岩相，测井数据为解释的密度、纵波速度和横波速度；
岩心数据和测井数据作为硬数据分配到划分的最邻近的网格节点上；
将工区范围内实际钻井上包含的信息称为硬数据或者称为条件数据，所述的硬数据包括了岩相、密度、纵波数据和横波速度；将分配了岩相的硬数据的网格称为已知网格或者称已知点；没有岩相数据的网格称为待估网格或又叫待估点；
在工区网格化和分配井数据后，所述的分配井数据为钻井数据和测井数据；
已知网格为分配有岩相数据的网格，
已知网格就是钻井穿过的网格，待估网格就是没有钻井穿过的网格；工区已知网格包含岩相、密度、纵波速度以及横波速度，待估网格没有任何数据；
③、赋予模拟工区初始属性值
根据统计的工区岩石弹性参数累积分布函数，对工区待估点赋予初始属 性值(指弹性参数)；
④、选择适当大小的数据样板
根据沉积相形态特征确定数据样板形状和大小；数据样板的形状可以根据沉积相的非均质性来定，非均质性强时可以采用二维的椭圆形或者三维的椭球状，非均质性较弱时可以采用二维长方形或者三维长方体；数据样板的大小没有硬性规定，一般可以采用5×7或者是5×5×7的数据样板大小；
⑤、反演
反演通过下面3个分步骤得到：
分步骤5.1：顺序访问模拟节点
先对工区所有网格随机产生伪随机数，该随机数的数值在0到1之间，然后以待估网格的数据样板大小为周围搜索条件点个数；每个待估点的数值等于伪随机数与条件点个数之和，对所有的待估点其进行从大到小排序，值较大的网格优先模拟，从而得到整个工区的模拟路径；
分步骤5.2：建议数据模式获取
如果进行初次反演，处理待估网格时，先要确定以待估点为中心的数据样板范围内的条件点的个数、相对待估点的位置以及岩相类型，形成以待估点为中心的以矢量形式呈现的数据事件；利用数据事件随机扫描步骤①中建立的训练图像，从中获得首个完全匹配的数据模式作为建议数据模式；建议数据模式表示一种沉积模式，是具有数据样板大小的岩相数据结构，数据样板中条件点的岩相数据结构在建议数据模式中可以完全体现；
如果进行迭代反演，处理待估网格时，先要根据待估点和待估点数据样板范围的条件点形成数据事件，利用数据事件随机扫描步骤①中建立的训练图像，从中获得若干完全匹配的数据模式作为候选数据模式；将候选数据模式按照中心点岩相类型排序，做成累积概率分布；
候选数据模式最多可包含从先到后训练图像中扫描到的前50个数据模式，当不到50个候选数据模式时，则用数据事件扫描完训练图像所有的数 据模式作为候选数据模式，根据候选数据模式中心岩相的不同计算不同岩相的比例，记做a 相1：a 相2：…：a 相m，再根据上一次弹性参数根据下述公式计算不同相的弹性参数的概率：
Figure PCTCN2020134088-appb-000010
公式中，b表示岩相弹性参数的概率，n表示待估点数据样板范围网格个数，ρ i表示待估点周围数据样板范围内第i个网格的密度值，
Figure PCTCN2020134088-appb-000011
表示待估点周围数据样板范围内第i个网格的纵波速度，
Figure PCTCN2020134088-appb-000012
表示待估点周围数据样板范围内第i个网格的横波速度，μ ρ表示该岩相密度均值，μ Vp表示该岩相纵波速度均值，μ Vs表示该岩相横波速度均值，
Figure PCTCN2020134088-appb-000013
表示该岩相密度方差，
Figure PCTCN2020134088-appb-000014
表示该岩相纵波速度方差，
Figure PCTCN2020134088-appb-000015
表示该岩相横波速度方差；
不同岩相的弹性参数的比例记做b 相1：b 相2：…：b 相m；
最后不同相的联合概率用(a 相1·b 相1)：(a 相2·b 相2):…:(a 相m·b 相m)计算，记做P 相1：P 相2：…：P 相m；将联合概率归一化处理，得到P’ 相1：P’ 相2：…：P’ 相m，将不同岩相的联合概率做成累积分布函数；
接着将候选数据模式按照中心网格的岩相进行排序，再随机抽样获取0到1之间的随机数，根据随机数在联合概率形成累积分布函数的位置和候选模式的累积概率分布，来确定建议数据模式；
分步骤5.3：建议弹性参数的获取
根据步骤①统计的不同岩石相弹性参数的累积分布函数，对建议数据模式中的各位置的岩相进行弹性参数的抽样；所述的弹性参数包括密度值和纵波速度值和横波速度值；
分步骤5.4：计算建议数据模式中建议弹性参数的合成记录
将建议弹性参数与预设入射角度基于下述Connolly公式计算反射系数；
Figure PCTCN2020134088-appb-000016
公式中，
Figure PCTCN2020134088-appb-000017
EI(θ)表示入射角为θ时的反射系数，θ代表预设入射角度，V p代表纵波速度，V s代表横波速度，ρ代表密度；
基于反射系数和地震地波，采用下述褶积公式计算初始正演模型：
g(m)＝w*EI(θ)
公式中，g(m)表示初始正演模型，m表示弹性参数，w表示地震子波；
分步骤5.5：选择和原始叠前数据匹配率最高的弹性参数
将上一次弹性参数和建议弹性参数通过分步骤5.3和5.4分别计算合成记录，前者记为g(m 上次)，后者记为g(m 建议)；这里上一次的弹性参数可以是初始弹性参数或可以是前一次迭代过程中的建议弹性参数；
采用如下公式将若干条合成地震道与实际地震记录比较，选择目标函数最小的地层模型参数和对应的建议数据模式作为最优反演结果：
Figure PCTCN2020134088-appb-000018
公式中，m代表选择的弹性参数，m 上次表示上一次的弹性参数，m 建议表示建议数据模型中建议弹性参数；
⑥、判断迭代终止
重复分步骤5.1-分步骤5.5，直到遍历所有网格，即所有待估点处模拟完成，就是实现了一次反演实现；
当迭代反演次数小于7次时，则重复步骤⑤，进行下一次迭代反演；当迭代反演第7次完成，则整个反演过程结束，得到多点地质统计学地震反演的模型，输出模拟结果。
在步骤②中；当一个网格中存在多个井数据时，选择离网格中心最近的数据点作为硬数据。
在分步骤5.1中；访问模拟节点的顺序通常是由井资料丰富区域逐步过渡到井资料较少区域，最后到无井区域。
在分步骤5.2中；岩石相条件数据包括井条件数据和已模拟节点的沉积相数据。
在分步骤5.2中；当不是初次迭代时，建议数据模式根据多点地质统计的先验地质信息和测井资料的弹性参数累积分布函数得到的联合概率共同决定。
在步骤③中；工区待估点赋予初始属性值可以是确定值，也可以是蒙特卡洛抽样随机获取。
在步骤①中；累积分布函数作为后期岩相约束下弹性参数抽样的对象。
在步骤③中；对工区待估点赋予初始属性值通过下面三个分步骤得到：
分步骤3.1：顺序访问模拟节点
访问模拟节点根据从一定顺序遍历所有网格，比如从下到上或从左到右或从前到后的原则；在网格中有岩相数据时，就不用赋予初始属性值；在网格中没有岩相数据时，则进行下一步；
分步骤3.2：根据蒙特卡洛抽样随机分配岩相；
根据步骤①中统计的不同岩相比例，采用蒙特卡洛抽样方法随机抽样，来对待估网格分配岩相，此时岩相只是临时保存在网格中，在下一分步骤根据岩相分配弹性参数后即可去掉；
分步骤3.3：根据蒙特卡洛抽样随机对不同岩相随机分配岩相对应的弹性参数；
根据分步骤3.2中网格的临时岩相数据，根据步骤①的不同岩相的不同弹性参数数据的累积概率分布，进行蒙特卡洛抽样获得该网格临时岩相的密度值、纵波速度值和横波速度值；
在赋予模拟工区初始属性值后，工区已知网格仍然是钻井穿过的网格，待估网格仍然是没有钻井穿过的网格；而岩相数据仅在已知网格上有，而弹性参数数据在工区所有网格中都有。
本公开实施例还包括如下具体实施例：
实施例1：参照图1为本实施例所述地震储层反演方法的实施步骤图，具体地：
(1)资料整理
根据工区河流相特征建立如图2所示的训练图像。建立工区砂泥岩岩石弹性参数的累积分布图。如图3、图4和图5所示，泥岩的密度范围是2.17-2.33g/cm3，纵波速度范围是3930-4272m/s，横波速度范围是2226-2496m/s；砂岩密度范围是2.27-2.42g/cm3，纵速度范围是4041-4528m/s，横速度范围是2300-2561m/s。砂泥岩的物性参数重叠交叉，各自又服从高斯分布。
(2)工区网格化和分配井数据
工区总长度1000m，总厚度45m。将模型网格划分为200×1×45，其中每个网格宽度5米，厚度1m。如图6所示，将实测井的相数据和岩石弹性参数分配到工区网格中。
(3)赋予模拟工区初始属性值
对工区待估点赋予初始属性值，该初始值用蒙特卡洛方法在统计的概率密度分布中抽取。
(4)选择适当大小的数据样板
根据工区河流相沉积相形态特征，选用5×1×5的数据样板
(5)确定访问待估网格的伪随机路径
先对工区所有网格随机生成伪随机数，然后数据样板大小范围搜索网格周围的条件点个数，两者数据之和就是待估点的值，对其进行从大到小排序，该顺序就是工区的模拟路径。以网格点(63，1，36)为例，数据样板范围内第一列为井数据，不参与随机数的赋值与排序，从左向右的第2、3两列整数为5表明每个网格周围数据样板范围包含5个井数据条件点，同样地，最右边的两列整数部分为0。
(6)获得待选数据模拟
以第一个网格点(63，1，36)为例，其中网格点(63，1，36)对应数据样板的中心点。用数据样板获取如图8所示的数据事件，开始扫描训练图像获得首个数据模式，该数据模式即建议数据模式(图9)。
(7)选择和原始合成记录匹配率最高的属性
选择主频为25HZ，总的时间长度为100ms，采样率为2ms的标准雷克子波作为地震子波。根据已建立的不同岩相的弹性参数累积分布图针对建议数据模式不同位置的岩相进行弹性参数的抽样，根据Connolly公式计算反射系数，将反射系数和地震子波褶积计算合成地震记录。将合成地震记录与原始地震记录比较并选择最优的一种模式与其对应的属性保留。对待估网格整体赋值后，转入下一个网格节点模拟直至所有网格节点被访问。图9、10、11、12和13分别为第一次迭代反演的相图、密度图、纵波速度图、横波速度图以及合成地震记录图。
(8)重新建立伪随机路径，开始后续的迭代反演。第二次迭代反演的相图、密度图、纵波速度图、横波速度图以及合成地震记录图如图14、15、16、17和18所示。与第一次迭代反演不同的是，在抽取前50个数据模式为数据模式后，根据中心点的相对数据模式排序，形成累积分布图，比如中心网格为泥岩相的数据模式有15个，为砂岩相的数据模式有35个。由于上一次迭代后网格密度值为2.32，纵波速度值为4320，横波速度值为2461，那么根据砂泥岩的概率密度函数计算，泥岩相和砂岩相之比为1:1，那么多点和地震的联合概率之比为3:7，先进行蒙特卡洛抽样判断中心网格是泥岩相还是砂岩相的数据模式，如抽样结果为0.4时，选择数据模式的累积分布图中对应第20个中心为砂岩相的数据模式作为建议数据模式。最后对建议数据模式不同位置的岩相进行弹性参数的抽样，根据Connolly公式计算反射系数，将反射系数和地震子波褶积计算合成地震记录。将合成地震记录与原始地震记录比较并选择最优的一种模式与其对应的属性保留。依次模拟所有网格。
(9)当迭代次数达到第七次时，第七次迭代反演的相图、密度图、纵波 速度图、横波速度图以及合成地震记录图如图19、20、21、22和23所示，输出反演结果，基于更新概率比率恒定理论的多点地质统计叠前反演完成。
上述未详细说明的部分均为现有技术。
以上所述仅为本公开实施例的较佳实施例，并不用以限制本公开实施例，凡在本公开实施例的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本公开实施例的保护范围之内。

Claims (14)

  1. 基于更新概率比率恒定理论的多点地质统计叠前反演方法，它包括如下具体步骤；
    ①、整理资料
    检查原始地震资料和井资料否齐全，其中地震资料为叠前地震数据，井资料包括钻井岩心和测井数据；根据井资料解释成果，可以确定不同井深对应的岩相以及不同岩相的弹性参数，弹性参数包括密度、纵波速度和横波速度；接着根据岩相数据建立不同岩相的相比例，根据弹性参数数据建立不同岩相弹性参数的累积分布函数；建立符合工区储层特征的训练图像；
    ②、工区网格化和分配井数据
    根据实际工区范围选择合适的网格尺寸大小，对工区进行网格划分，建立网格模型；根据各井的平面位置和井深将岩心数据和测井数据网格化；岩心数据即为岩相，测井数据为解释的密度、纵波速度和横波速度；
    岩心数据和测井数据作为硬数据分配到划分的最邻近的网格节点上；
    将工区范围内实际钻井上包含的信息称为硬数据或者称为条件数据，所述的硬数据包括了岩相、密度、纵波数据和横波速度；将分配了岩相的硬数据的网格称为已知网格或者称已知点；没有岩相数据的网格称为待估网格或又叫待估点；
    在工区网格化和分配井数据后，所述的分配井数据为钻井数据和测井数据；
    已知网格为分配有岩相数据的网格，
    已知网格就是钻井穿过的网格，待估网格就是没有钻井穿过的网格；工区已知网格包含岩相、密度、纵波速度以及横波速度，待估网格没有任何数据；
    ③、赋予模拟工区初始属性值
    根据统计的工区岩石弹性参数累积分布函数，对工区待估点赋予初始属性值，所述初始属性值指弹性参数；
    ④、选择适当大小的数据样板
    根据沉积相形态特征确定数据样板形状和大小；数据样板的形状可以根据沉积相的非均质性来定，非均质性强时可以采用二维的椭圆形或者三维的椭球状，非均质性较弱时可以采用二维长方形或者三维长方体；数据样板的大小没有硬性规定，一般可以采用5×7或者是5×5×7的数据样板大小；
    ⑤、反演
    反演过程包括：
    分步骤5.1：顺序访问模拟节点；
    分步骤5.2：建议数据模式获取；
    分步骤5.3：建议弹性参数的获取；
    分步骤5.4：计算建议数据模式中建议弹性参数的合成记录；
    分步骤5.5：选择和原始叠前数据匹配率最高的弹性参数；
    ⑥、判断迭代终止后输出模拟结果。
  2. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在步骤①中；累积分布函数作为后期岩相约束下弹性参数抽样的对象。
  3. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在步骤②中；当一个网格中存在多个井数据时，选择离网格中心最近的数据点作为硬数据。
  4. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在步骤③中；工区待估点赋予初始属性值可以是确定值，也可以是蒙特卡洛抽样随机获取。
  5. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在步骤③中；
    对工区待估点赋予初始属性值通过下面三个分步骤得到：
    分步骤3.1：顺序访问模拟节点
    访问模拟节点根据从一定顺序遍历所有网格，所述一定顺序为从下到上或从左到右或从前到后的原则；在网格中有岩相数据时，就不用赋予初始属性值；在网格中没有岩相数据时，则进行下一步；
    分步骤3.2：根据蒙特卡洛抽样随机分配岩相；
    根据步骤①中统计的不同岩相比例，采用蒙特卡洛抽样方法随机抽样，来对待估网格分配岩相，此时岩相只是临时保存在网格中，在下一分步骤根据岩相分配弹性参数后即可去掉；
    分步骤3.3：根据蒙特卡洛抽样随机对不同岩相随机分配岩相对应的弹性参数；
    根据分步骤3.2中网格的临时岩相数据，根据步骤①的不同岩相的不同弹性参数数据的累积概率分布，进行蒙特卡洛抽样获得该网格临时岩相的密度值、纵波速度值和横波速度值；
    在赋予模拟工区初始属性值后，工区已知网格仍然是钻井穿过的网格，待估网格仍然是没有钻井穿过的网格；而岩相数据仅在已知网格上有，而弹性参数数据在工区所有网格中都有。
  6. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，所述分步骤5.1包括：
    先对工区所有网格随机产生伪随机数，该随机数的数值在0到1之间，然后以待估网格的数据样板大小为周围搜索条件点个数；每个待估点的数值等于伪随机数与条件点个数之和，对所有的待估点其进行从大到小排序，值较大的网格优先模拟，从而得到整个工区的模拟路径。
  7. 根据权利要求6所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在分步骤5.1中；访问模拟节点的顺序通常是由井资料丰富区域逐步过渡到井资料较少区域，最后到无井区域。
  8. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，所述分步骤5.2包括：
    如果进行初次反演，处理待估网格时，先要确定以待估点为中心的数据样板范围内的条件点的个数、相对待估点的位置以及岩相类型，形成以待估点为中心的以矢量形式呈现的数据事件；利用数据事件随机扫描步骤①中建立的训练图像，从中获得首个完全匹配的数据模式作为建议数据模式；建议数据模式表示一种沉积模式，是具有数据样板大小的岩相数据结构，数据样板中条件点的岩相数据结构在建议数据模式中可以完全体现；
    如果进行迭代反演，处理待估网格时，先要根据待估点和待估点数据样板范围的条件点形成数据事件，利用数据事件随机扫描步骤①中建立的训练图像，从中获得若干完全匹配的数据模式作为候选数据模式；将候选数据模式按照中心点岩相类型排序，做成累积概率分布；
    候选数据模式最多可包含从先到后训练图像中扫描到的前50个数据模式，当不到50个候选数据模式时，则用数据事件扫描完训练图像所有的数据模式作为候选数据模式，根据候选数据模式中心岩相的不同计算不同岩相的比例，记做a 相1：a 相2：…：a 相m，再根据上一次弹性参数根据下述公式计算不同相的弹性参数的概率：
    Figure PCTCN2020134088-appb-100001
    公式中，b表示岩相弹性参数的概率，n表示待估点数据样板范围网格个数，ρ i表示待估点周围数据样板范围内第i个网格的密度值，
    Figure PCTCN2020134088-appb-100002
    表示待估点周围数据样板范围内第i个网格的纵波速度，
    Figure PCTCN2020134088-appb-100003
    表示待估点周围数据样板范围内第i个网格的横波速度，μ ρ表示该岩相密度均值，μ Vp表示该岩相纵波速度均值，μ Vs表示该岩相横波速度均值，
    Figure PCTCN2020134088-appb-100004
    表示该岩相密度方差，
    Figure PCTCN2020134088-appb-100005
    表示该岩相纵波速度方差，
    Figure PCTCN2020134088-appb-100006
    表示该岩相横波速度方差；
    不同岩相的弹性参数的比例记做b 相1：b 相2：…：b 相m；
    最后不同相的联合概率用(a 相1·b 相1)：(a 相2·b 相2):…:(a 相m·b 相m)计算，记做P 相1：P 相2：…：P 相m；将联合概率归一化处理，得到P’ 相1：P’ 相2：…：P’ 相m，将不同岩相的联合概率做成累积分布函数；
    接着将候选数据模式按照中心网格的岩相进行排序，再随机抽样获取0到1之间的随机数，根据随机数在联合概率形成累积分布函数的位置和候选模式的累积概率分布，来确定建议数据模式。
  9. 根据权利要求8所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在分步骤5.2中；岩石相条件数据包括井条件数据和已模拟节点的沉积相数据。
  10. 根据权利要求8所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，在分步骤5.2中；当不是初次迭代时，建议数据模式根据多点地质统计的先验地质信息和测井资料的弹性参数累积分布函数得到的联合概率共同决定。
  11. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，所述分步骤5.3包括：
    根据步骤①统计的不同岩石相弹性参数的累积分布函数，对建议数据模式中的各位置的岩相进行弹性参数的抽样；所述的弹性参数包括密度值和纵波速度值和横波速度值。
  12. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，所述分步骤5.4包括：
    将建议弹性参数与预设入射角度基于下述Connolly公式计算反射系数；
    Figure PCTCN2020134088-appb-100007
    公式中，
    Figure PCTCN2020134088-appb-100008
    EI(θ)表示入射角为θ时的反射系数，θ代表预设入射角度，V p代表纵波速度，V s代表横波速度，ρ代表密度；
    基于反射系数和地震地波，采用下述褶积公式计算初始正演模型：
    g(m)＝w*EI(θ)
    公式中，g(m)表示初始正演模型，m表示弹性参数，w表示地震子波。
  13. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，所述分步骤5.5包括：
    将上一次弹性参数和建议弹性参数通过分步骤5.3和5.4分别计算合成记录，前者记为g(m上次)，后者记为g(m建议)；这里上一次的弹性参数可以是初始弹性参数或可以是前一次迭代过程中的建议弹性参数；
    采用如下公式将若干条合成地震道与实际地震记录比较，选择目标函数最小的地层模型参数和对应的建议数据模式作为最优反演结果：
    Figure PCTCN2020134088-appb-100009
    公式中，m代表选择的弹性参数，m 上次表示上一次的弹性参数，m 建议表示建议数据模型中建议弹性参数。
  14. 根据权利要求1所述的基于更新概率比率恒定理论的多点地质统计叠前反演方法，其中，⑥、判断迭代终止的过程包括：
    重复分步骤5.1-分步骤5.5，直到遍历所有网格，即所有待估点处模拟完成，就是实现了一次反演实现；
    当迭代反演次数小于7次时，则重复步骤⑤，进行下一次迭代反演；当迭代反演第7次完成，则整个反演过程结束，得到多点地质统计学地震反演的模型，输出模拟结果。
PCT/CN2020/134088 2020-01-21 2020-12-04 基于更新概率比率恒定理论的多点地质统计叠前反演方法 WO2021147529A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US17/315,503 US11131784B2 (en) 2020-01-21 2021-05-10 Multi-point geostatistical prestack inversion method based on renewal probability ratio constant theory

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202010068937.8A CN111273348B (zh) 2020-01-21 2020-01-21 基于更新概率比率恒定理论的多点地质统计叠前反演方法
CN202010068937.8 2020-01-21

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US17/315,503 Continuation US11131784B2 (en) 2020-01-21 2021-05-10 Multi-point geostatistical prestack inversion method based on renewal probability ratio constant theory

Publications (1)

Publication Number Publication Date
WO2021147529A1 true WO2021147529A1 (zh) 2021-07-29

Family

ID=70997628

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2020/134088 WO2021147529A1 (zh) 2020-01-21 2020-12-04 基于更新概率比率恒定理论的多点地质统计叠前反演方法

Country Status (3)

Country Link
US (1) US11131784B2 (zh)
CN (1) CN111273348B (zh)
WO (1) WO2021147529A1 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114329974A (zh) * 2021-12-29 2022-04-12 武汉大学 基于蒙特卡洛模拟的城市供水管网地震损伤评估方法
CN114779335A (zh) * 2022-03-31 2022-07-22 吉林大学 基于各向异性全变分约束的弹性波直接包络反演方法
CN117388919A (zh) * 2023-08-29 2024-01-12 河海大学 一种基于dnn预测致密油储层横波速度的方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111273348B (zh) * 2020-01-21 2021-02-05 长江大学 基于更新概率比率恒定理论的多点地质统计叠前反演方法
CN113406702B (zh) * 2021-05-10 2023-04-07 中国石油天然气集团有限公司 基于相约束可变网格的地质统计学反演方法及装置
CN114608527B (zh) * 2022-03-08 2024-03-26 内蒙古银宏能源开发有限公司 一种基于煤矿开采地表沉陷特征的D-InSAR和UAV摄影测量融合方法
CN114779324B (zh) * 2022-03-11 2023-07-25 西南交通大学 一种基于深度学习的瑞雷波频散曲线反演方法
CN115421181B (zh) * 2022-07-27 2023-10-20 北京超维创想信息技术有限公司 一种基于深度学习的三维地质模型相控属性建模方法
CN115903034B (zh) * 2022-10-28 2023-08-08 中国海洋大学 基于弹性参数概率分布的水合物与游离气共存层识别方法
CN116976751B (zh) * 2023-08-28 2024-02-23 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) 一种基于时空数据融合的煤矿采空区勘察系统
CN117075212B (zh) * 2023-10-16 2024-01-26 吉林大学 一种隧道磁共振裂隙结构成像方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100157730A1 (en) * 2008-12-23 2010-06-24 Schlumberger Technology Corporation Method of subsurface imaging using microseismic data
CN104850682A (zh) * 2015-04-17 2015-08-19 长江大学 基于位置的多点地质统计学建模方法
CN106154323A (zh) * 2015-04-01 2016-11-23 中国石油化工股份有限公司 基于地震拓频处理的相控随机反演薄储层预测方法
CN108645994A (zh) * 2018-04-25 2018-10-12 中国石油大学(北京) 一种基于多点地质统计学的地质随机反演方法及装置
CN108663710A (zh) * 2017-03-30 2018-10-16 中国石油化工股份有限公司 宽方位地震数据处理一体化成像反演方法及系统
CN108931811A (zh) * 2018-05-17 2018-12-04 长江大学 基于多点地质统计的地震储层反演方法
CN110031896A (zh) * 2019-04-08 2019-07-19 中国石油天然气集团有限公司 基于多点地质统计学先验信息的地震随机反演方法及装置
CN110031895A (zh) * 2019-03-11 2019-07-19 西安科技大学 一种基于图像缝合的多点地质统计学随机反演方法及装置
CN111273348A (zh) * 2020-01-21 2020-06-12 长江大学 基于更新概率比率恒定理论的多点地质统计叠前反演方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7254091B1 (en) * 2006-06-08 2007-08-07 Bhp Billiton Innovation Pty Ltd. Method for estimating and/or reducing uncertainty in reservoir models of potential petroleum reservoirs
US8729903B2 (en) * 2009-11-09 2014-05-20 Exxonmobil Upstream Research Company Method for remote identification and characterization of hydrocarbon source rocks using seismic and electromagnetic geophysical data
WO2012015542A1 (en) * 2010-07-27 2012-02-02 Exxonmobil Upstream Research Company Inverting geophysical data for geological parameters or lithology
EP2856373A4 (en) * 2012-05-24 2016-06-29 Exxonmobil Upstream Res Co SYSTEM AND METHOD FOR PREDICTING THE RESISTANCE OF A ROCK

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100157730A1 (en) * 2008-12-23 2010-06-24 Schlumberger Technology Corporation Method of subsurface imaging using microseismic data
CN106154323A (zh) * 2015-04-01 2016-11-23 中国石油化工股份有限公司 基于地震拓频处理的相控随机反演薄储层预测方法
CN104850682A (zh) * 2015-04-17 2015-08-19 长江大学 基于位置的多点地质统计学建模方法
CN108663710A (zh) * 2017-03-30 2018-10-16 中国石油化工股份有限公司 宽方位地震数据处理一体化成像反演方法及系统
CN108645994A (zh) * 2018-04-25 2018-10-12 中国石油大学(北京) 一种基于多点地质统计学的地质随机反演方法及装置
CN108931811A (zh) * 2018-05-17 2018-12-04 长江大学 基于多点地质统计的地震储层反演方法
CN110031895A (zh) * 2019-03-11 2019-07-19 西安科技大学 一种基于图像缝合的多点地质统计学随机反演方法及装置
CN110031896A (zh) * 2019-04-08 2019-07-19 中国石油天然气集团有限公司 基于多点地质统计学先验信息的地震随机反演方法及装置
CN111273348A (zh) * 2020-01-21 2020-06-12 长江大学 基于更新概率比率恒定理论的多点地质统计叠前反演方法

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114329974A (zh) * 2021-12-29 2022-04-12 武汉大学 基于蒙特卡洛模拟的城市供水管网地震损伤评估方法
CN114779335A (zh) * 2022-03-31 2022-07-22 吉林大学 基于各向异性全变分约束的弹性波直接包络反演方法
CN117388919A (zh) * 2023-08-29 2024-01-12 河海大学 一种基于dnn预测致密油储层横波速度的方法
CN117388919B (zh) * 2023-08-29 2024-05-28 河海大学 一种基于dnn预测致密油储层横波速度的方法

Also Published As

Publication number Publication date
US11131784B2 (en) 2021-09-28
US20210263176A1 (en) 2021-08-26
CN111273348A (zh) 2020-06-12
CN111273348B (zh) 2021-02-05

Similar Documents

Publication Publication Date Title
WO2021147529A1 (zh) 基于更新概率比率恒定理论的多点地质统计叠前反演方法
CN110031896B (zh) 基于多点地质统计学先验信息的地震随机反演方法及装置
CN113759424B (zh) 基于频谱分解和机器学习的岩溶储层充填分析方法和系统
WO2018129844A1 (zh) 一种地震绕射波分离方法和装置
WO2020123084A1 (en) Machine learning-augmented geophysical inversion
US8666149B2 (en) Method for editing a multi-point facies simulation
CN104755960B (zh) 基于盆地建模对用于处理地震数据的速度模型的改进
WO2017007924A1 (en) Improved geobody continuity in geological models based on multiple point statistics
EP1782328A2 (en) Multiple-point statistic (mps) simulation with enhanced computational efficiency
CN111505713B (zh) 基于多点地质统计的叠前地震反演方法
CN108931811A (zh) 基于多点地质统计的地震储层反演方法
CN110927793B (zh) 一种基于序贯随机模糊模拟的储层预测方法及系统
US6324478B1 (en) Second-and higher-order traveltimes for seismic imaging
US20210405250A1 (en) Conditioning of Surface-Based Geologic Models
CN115407407A (zh) 针对碳酸盐岩古溶洞及其充填的三维地质模型构建方法
Strebelle Sequential simulation for modeling geological structures from training images
CN108537883B (zh) 一种基于MapReduce框架的地质建模方法
CN115880455A (zh) 基于深度学习的三维智能插值方法
CN114153002A (zh) 储层天然裂缝三维地质建模方法、装置、电子设备及介质
Huang et al. A Novel Method of 3D Multipoint Geostatistical Inversion Using 2D Training Images
EP3320450B1 (en) Improved geobody continuity in geological models based on multiple point statistics
CN115469361B (zh) 一种碎屑岩地层三维地质建模方法
CN117407841B (zh) 一种基于优化集成算法的页岩层理缝预测方法
WO2024087800A1 (zh) 基于常规测井资料的高频旋回碳酸盐岩识别方法
CN110488350B (zh) 基于卷积神经网络的地震反演大数据生成方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 20915138

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20915138

Country of ref document: EP

Kind code of ref document: A1

122 Ep: pct application non-entry in european phase

Ref document number: 20915138

Country of ref document: EP

Kind code of ref document: A1