A numerical simulation method for seismic response of organic-rich shale
Through the image processing and the establishment of a lattice spring model of organic mud rock heart samples, the existing models fail to fully consider the problem of organic matter and clay mineral systems, improve the accuracy of earthquake response modeling, and develop for oil and gas exploration and development Provide important theoretical support.
Patent Information
- Application Number
- CN202210862994.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2022-07-21
- Publication Date
- 2025-01-28
- Estimated Expiration
- 2042-07-21
AI Technical Summary
The existing seismic petrophysical models fail to fully consider organic matter and clay minerals as an interrelated system, resulting in insufficient research on seismic responses to organic matter-rich mud shale.
Through information collection and image processing of organic mud rock heart samples, the rock physical model is generated, and a lattice spring model is established to simulate the earthquake response and fully consider the micro -structure of organic clay composite.
It improves the accuracy of the wealthy mud rock earthquake response, which can more accurately describe the quantitative relationship between the micro -structure of the rock and the earthquake response, and provides an important theoretical basis for the exploration and development of land shale oil.
Smart Images

Figure CN115267892B_ABST
Abstract
Description
Technical Field
[0001] The invention provides a numerical simulation method for seismic response of organic-rich shale, belonging to the field of oil and gas geophysics. Background Art
[0002] The organic-rich shale in continental lake basins is different from the marine shale in North America, which has a high degree of exploration and development. It has the characteristics of rapid shale phase change, complex structure, low formation energy, strong reservoir heterogeneity, heavy oil, high viscosity and poor fluidity, which leads to huge challenges in the exploration and development of continental shale oil. Therefore, it is imperative to strengthen relevant basic research.
[0003] Seismic response characteristics are one of the key research contents. This study is a direct bridge to obtain reservoir physical parameters (organic matter, clay content, porosity, etc.) and mechanical characteristics (elastic modulus, brittleness index, etc.) through seismic attributes. As the first step, it is necessary to build a seismic rock physics model and establish a quantitative relationship between rock microstructure parameters and macroscopic seismic response. The seismic response is controlled by its own structural characteristics. How to fully consider the spatial distribution characteristics of the main constituent minerals, organic matter, pore structure, etc., and establish a set of theoretical models that more accurately describe the relevant structural characteristics from a microscopic scale is the key to studying the seismic rock physics characteristics and response mechanism of continental organic-rich shale.
[0004] Organic matter and clay minerals are two important components of mudstone. The organic clay complex formed by the two during the sedimentary evolution process is the natural parent material for oil and gas generation, and plays an important role in the enrichment and preservation of organic matter in mudstone. Existing seismic rock physics models mostly focus on the characterization of organic matter and pore structure, and the research on clay minerals and the relationship between clay minerals and organic matter is still insufficient, and the importance of studying organic matter and clay minerals as an interrelated system has not been fully considered. Summary of the invention
[0005] The purpose of the present invention is to provide a numerical simulation method for the simulation modeling study of the seismic response of organic-rich shale on the basis of fully considering the microstructure of the organic clay complex.
[0006] To achieve the above object, the present invention provides a numerical simulation method for seismic response of organic-rich shale, comprising the following steps:
[0007] 1) acquiring a digital two-dimensional image of the core sample and capturing the main mineral components and spatial distribution characteristics of the core sample by performing information collection and image processing on the organic-rich shale core sample;
[0008] 2) generating a rock physics model of the organic-rich shale sample according to the main mineral components and spatial distribution characteristics of the core sample obtained in step 1;
[0009] 3) According to the rock physics model generated in step 2, a corresponding lattice spring model is established for seismic response simulation modeling.
[0010] The step 1) comprises:
[0011] 11) Collecting information and processing images of the core sample includes converting the digital core two-dimensional image of the core sample into a grayscale image using the following color-to-grayscale value formula (1):
[0012] I Gray =I Red 0.299+I Green 0.587+I Blue 0.114 (2)
[0013] In the formula, I Gray Represents gray value; I Red ,I Green ,I Blue Represent the red, green, and blue component values respectively;
[0014] 12) Determine the main minerals and spatial distribution characteristics of the core samples based on the differences in X-ray transmittance of different media;
[0015] 13) Reconstructing the digital core two-dimensional image of the core sample into a three-dimensional image by combining the grayscale image obtained in step 11) and the main mineral and spatial distribution characteristics of the core sample obtained in step 12), and naming the three-dimensional image digital core.
[0016] In the step 12), since pyrite and pore structure have a higher difference in absorption coefficient, while organic matter, brittle minerals, clay minerals, organic clay complexes, etc. have similar absorption coefficients, the medium of the core sample is preliminarily divided into three groups of media: irregular pore structure, pyrite and preliminary rock skeleton.
[0017] The information collection method of the core sample in step 1) includes high-precision CT imaging, comprehensive mineral analysis electron microscopy and / or scanning electron microscopy.
[0018] The step 2) comprises:
[0019] 21) Counting the distribution characteristics of organic matter, organoclay complexes, and clay minerals in the preliminary rock skeleton in step 12, determining a random medium statistical model, and generating a three-dimensional random medium image based on the digital core size;
[0020] 22) performing Boolean operation on the three-dimensional random medium image, performing binarization image processing on the three-dimensional random medium image after the Boolean operation, and performing watershed transformation, and using the three-dimensional random medium image to simulate and divide various structural materials of the preliminary rock skeleton;
[0021] 23) Estimate the elastic material parameters of each structural material obtained in step 22).
[0022] The structural materials described in step 22) include organic matter, viscosity minerals, agglomerates of organic clay complexes and the final rock skeleton.
[0023] The step 2) is specifically as follows:
[0024] 21) Statistical analysis of the distribution characteristics of organic matter, organic clay complexes, and clay minerals in core samples, determine a reasonable random medium statistical model, and generate a three-dimensional random medium image based on the size of the digital core sample, denoted as V sample , using a mixed random medium three-dimensional model:
[0025]
[0026] Where a, b, and c are the autocorrelation lengths of the medium in the x, y, and z directions, and r is called the roughness factor. These parameters are estimated by inverting the scanned image;
[0027] 22) Perform Boolean operations on three-dimensional random medium images,
[0028] V rest =V sample -(V pore ∪V pyrite ) (4)
[0029] For the three-dimensional image V after Boolean operation rest Binarized image processing and watershed transformation were performed to divide the three-dimensional image into clusters representing organic matter, viscosity minerals, organic clay complexes, and the final rock skeleton, and the corresponding geometric information vector files were obtained, which were recorded as V1, ..., V i ,……,V N , V matrix , where the subscript i represents the i-th cluster, N represents the number of identified clusters, and V matrix For the final rock skeleton;
[0030] 23) Calculate the irregular pore structure V that makes up the digital core sample pore , Pyrite V pyrite , the clusters generated based on random media are V1, ..., V i ,……,V N , and finally the rock skeleton Vmatrix Elastic material parameters,
[0031] These parameters satisfy the following formula:
[0032] V sample =V1+…+V i +…+V N +V matrix +V pore +V pyrite (5)
[0033] Assuming that each group of media satisfies the isotropic characteristics, each group only needs two independent elastic material parameters, where the bulk modulus K and the shear modulus μ are selected, where the pore structure V pore With pyrite V pyrite The elastic material parameters of the rock are considered to be composed of one component; while the agglomerates and the final rock skeleton may be composed of multiple components, although the final rock skeleton V matrix It is composed of a variety of minerals including clay minerals and brittle minerals, but here we can only focus on the overall material characteristics of the final rock skeleton, which is also described by two independent material elastic parameters; assuming that the components of the agglomerate are mainly composed of clay and organic matter, the proportions of the two in the i-th agglomerate are recorded as and Both meet
[0034]
[0035] When hour, Indicates that the mass is organic matter. hour, Indicates that the mass is pure clay. and When , it means that the agglomerate is an organic clay complex, and each agglomerate has and The elastic material parameters of each agglomerate are estimated based on the average value of the Hashin-Shtrikman boundary theory according to the clay and organic matter content of the sample. The volume model and shear modulus of the i-th agglomerate are expressed as K i With μ i express:
[0036]
[0037]
[0038]
[0039]
[0040]
[0041]
[0042] The step 3) comprises:
[0043] 31) Generate randomly and closely arranged lattice nodes according to the size of the digital core, determine the lattice nodes representing each group of minerals, and classify and number them;
[0044] 32) generating a lattice spring network based on a nearest neighbor search algorithm, and determining input parameters of each lattice node in the lattice spring model;
[0045] 33) Based on the ultrasonic experimental data, the micro spring stiffness parameters representing the final rock skeleton are identified, and the numerical simulation of seismic response is carried out based on the lattice spring model, that is, the lattice spring model for numerical simulation of seismic response of organic-rich shale taking into account the microstructure of the organic clay complex is obtained.
[0046] The step 3) is specifically as follows:
[0047] 31) According to the digital core size V sample Generate randomly closely arranged lattice nodes and determine the lattice nodes representing each group of minerals. A total of N+3 three-dimensional geometric information files (i.e., .stl files) are generated by steps 1) and 2). The respective geometric spaces to which the lattice nodes belong are determined in a cyclic manner, and the irregular pore structure V is deleted. pore The lattice nodes in the represented geometric space are classified and numbered;
[0048] 32) Generate a lattice spring network based on the nearest neighbor search algorithm and determine the representative blocks (V1, ..., V i ,……,V N ) and pyrite V pyrite The input parameters of the lattice node. The elastic parameters of the lattice node mainly include the normal spring stiffness k n and the tangential spring stiffness k s , the micro spring parameters and macro material properties satisfy the following relationship:
[0049]
[0050]
[0051]
[0052] Where α is the microstructure coefficient, l j is the length of the jth lattice spring, V mis the volume of each three-dimensional geometric space, K and μ are the bulk modulus and shear modulus of the macroscopic material represented by each three-dimensional geometric space, and the macroscopic elastic material parameters of the mass are obtained in step 2) formulas (6)-(10), pyrite V pyrite The macro elastic material parameters of the rock can be obtained by consulting the rock mineral manual. matrix The macroscopic elastic material parameters are to be determined;
[0053] 33) Perform ultrasonic testing on the core sample in step 1), and identify the V representing the final rock skeleton by integrating the ultrasonic test data. matrix The micro spring stiffness parameter can be obtained through ultrasonic experiments. The longitudinal wave v of ultrasonic waves propagating in the medium P and the shear wave velocity v s , and then the macroscopic bulk modulus K of the sample can be estimated by the following formula: sample and shear modulus μ sample
[0054]
[0055]
[0056] The objective function of constructing numerical simulation data and ultrasonic experimental data based on the least squares method is shown in the following formula:
[0057]
[0058] Where s represents the input parameter data set to be calibrated, including the final rock skeleton V matrix The micro spring stiffness parameter, that is, the normal spring stiffness k n and the tangential spring stiffness k s , K(s) and μ(s) represent the simulation results obtained with different s. After multiple iterations of calibration, the optimal input parameter data set s is obtained, that is, the lattice spring model for numerical simulation of seismic response of organic-rich shale taking into account the microstructure of the organic clay complex is obtained.
[0059] Wherein, the three-dimensional geometric information file is an stl file.
[0060] Among them, the numerical simulation method for the seismic response of organic-rich shale is a numerical simulation method for the seismic response of organic-rich shale that can fully consider the microstructure of the organic clay complex and take into account the microstructure of the organic clay complex.
[0061] Due to the adoption of the above technical scheme, the present invention has the following advantages: the present invention integrates advanced image processing technology, random medium theory and discrete modeling methods, and proposes a numerical simulation method for the seismic response of organic-rich shale taking into account the microstructure of organic-clay complex. The model fully considers the microstructural characteristics of organic-clay complex in organic-rich shale, and innovatively adopts lattice spring model for seismic response simulation of such complex medium, which can effectively improve the modeling accuracy and fully consider the abnormal wave field characteristics induced by small-scale microstructure, so as to provide the necessary theoretical basis for the identification and prediction of oil and gas reservoir enrichment areas based on elastic wave theory. This research will have very important scientific significance for the exploration and development of continental shale oil in my country. BRIEF DESCRIPTION OF THE DRAWINGS
[0062] Figure 1 It is a schematic diagram of the process of the present invention. DETAILED DESCRIPTION
[0063] The present invention is further described in detail below in conjunction with specific embodiments, and the examples provided are only for illustrating the present invention, rather than for limiting the scope of the present invention. The examples provided below can be used as a guide for further improvements by those of ordinary skill in the art, and do not constitute a limitation of the present invention in any way.
[0064] The experimental methods in the following examples, unless otherwise specified, are all conventional methods, and are performed according to the techniques or conditions described in the literature in the field or according to the product instructions. The materials, reagents, etc. used in the following examples, unless otherwise specified, can all be obtained from commercial channels.
[0065] The present invention adopts the following technical scheme: a numerical simulation method for seismic response of organic-rich shale taking into account the microstructure of organic clay complex, comprising the following steps:
[0066] 1) Based on high-precision CT imaging and image processing technology, the spatial distribution characteristics of the main mineral components and pore structure of the core samples are "truly" captured;
[0067] 2) Generate rock physics models of organic-rich shale samples that take into account the microstructure of organoclay complexes;
[0068] 3) Based on the rock physics model, a corresponding lattice spring model is established for seismic response simulation modeling.
[0069] The numerical simulation method for seismic response of organic-rich shale taking into account the microstructure of the organoclay complex, preferably, the specific process of step 1) is as follows:
[0070] ① The digital core two-dimensional image is converted into a grayscale image using the following color-to-grayscale value formula (1):
[0071] I Gray =I Red 0.299+I Green 0.587+I Blue 0.114 (19)
[0072] In the formula, I Gray Represents gray value; I Red ,I Green ,I Blue Represent the red, green, and blue component values respectively;
[0073] ② According to the difference of X-ray transmittance of different media, the spatial distribution characteristics of the main minerals in the core samples are determined. Since pyrite and pore structure have high absorption coefficient differences, and organic matter, brittle minerals, clay minerals, organic clay complexes, etc. have similar absorption coefficients, the digital cores are initially divided into three groups of media: irregular pore structure, pyrite, and preliminary rock skeleton (the preliminary division result of the rock skeleton is distinguished from the final rock skeleton generated in step 2).
[0074] ③ The Marching Cubes algorithm is used to reconstruct the two-dimensional images of irregular pore structure and pyrite into three-dimensional images, and the corresponding geometric information vector files (generally in .stl format) are obtained, which are recorded as V pore With V pyrite .
[0075] The numerical simulation method for seismic response of organic-rich shale taking into account the microstructure of the organoclay complex, preferably, the specific process of step 2) is as follows:
[0076] ① Comprehensively analyze the results of mineral analysis electron microscope, scanning electron microscope and other microscopic imaging results, statistically analyze the distribution characteristics of organic matter, organic clay complex and clay minerals in core samples, determine a reasonable random medium statistical model, and generate a three-dimensional random medium image based on the size of the digital core sample, denoted as V sample , using a mixed random medium three-dimensional model:
[0077]
[0078] Where a, b, and c are the autocorrelation lengths of the medium in the x, y, and z directions, and r is called the roughness factor. These parameters are estimated by inverting the scanned image.
[0079] ② Perform Boolean operations on three-dimensional random medium images.
[0080] V rest =V sample -(V pore ∪V pyrite) (twenty one)
[0081] For the three-dimensional image V after Boolean operation rest Binarized image processing and watershed transformation are performed to divide the three-dimensional image into clumps representing organic matter, viscosity minerals, organic clay complexes, etc., and the final rock skeleton, and the corresponding geometric information vector files (generally in .stl format) are obtained, which are recorded as V1, ..., V i ,……,V N , V matrix , where the subscript i represents the i-th cluster, N represents the number of identified clusters, and V matrix For the final final rock skeleton;
[0082] ③ Calculate the irregular pore structure V that makes up the digital core sample pore , Pyrite V pyrite , the clusters generated based on random media are V1, ..., V i ,……,V N , and finally the rock skeleton V matrix Elastic material parameters, which satisfy the following formula:
[0083] V sample =V1+...+V i +...+V N +V matrix +V pore +V pyrite (twenty two)
[0084] Assuming that each group of media satisfies the isotropic characteristics, each group only needs two independent elastic material parameters, where the bulk modulus K and the shear modulus μ are selected, where the pore structure V pore With pyrite V pyrite The elastic material parameters of the rock are considered to be composed of one component; while the agglomerates and the final rock skeleton may be composed of multiple components, although the final rock skeleton V matrix It is composed of a variety of minerals including clay minerals and brittle minerals, but here we can only focus on the overall material characteristics of the final rock skeleton, which is also described by two independent material elastic parameters; assuming that the components of the agglomerate are mainly composed of clay and organic matter, the proportions of the two in the i-th agglomerate are recorded as and Both meet
[0085]
[0086] When hour, Indicates that the mass is organic matter. hour, Indicates that the mass is pure clay. and When , it means that the agglomerate is an organic clay complex, and each agglomerate has and The elastic material parameters of each agglomerate are estimated based on the average value of the Hashin-Shtrikman boundary theory according to the clay and organic matter content of the sample. The volume model and shear modulus of the i-th agglomerate are expressed as K i With μ i express:
[0087]
[0088]
[0089]
[0090]
[0091]
[0092]
[0093] The numerical simulation method for seismic response of organic-rich shale taking into account the microstructure of the organoclay complex, preferably, the specific process of step 3) is as follows:
[0094] ①According to the digital core size V sample Generate randomly closely arranged lattice nodes and determine the lattice nodes representing each group of minerals. A total of N+3 three-dimensional geometric information files (i.e., .stl files) are generated by steps 1) and 2). The respective geometric spaces to which the lattice nodes belong are determined in a cyclic manner, and the irregular pore structure V is deleted. pore The lattice nodes in the represented geometric space are classified and numbered;
[0095] ② Generate a lattice spring network based on the nearest neighbor search algorithm and determine the representative blocks (V1, ..., V i ,……,V N ) and pyrite V pyrite The input parameters of the lattice node. The elastic parameters of the lattice node mainly include the normal spring stiffness k n and the tangential spring stiffness k s , the micro spring parameters and macro material properties satisfy the following relationship:
[0096]
[0097]
[0098]
[0099] Where α is the microstructure coefficient, l j is the length of the jth lattice spring, V m is the volume of each three-dimensional geometric space, K and μ are the bulk modulus and shear modulus of the macroscopic material represented by each three-dimensional geometric space, and the macroscopic elastic material parameters of the mass are obtained in step 2) formulas (6)-(10), pyrite V pyrite The macro elastic material parameters of the rock can be obtained by consulting the rock mineral manual. matrix The macroscopic elastic material parameters are to be determined;
[0100] ③Combining ultrasonic experimental data, identify the final rock skeleton V matrix The micro spring stiffness parameter can be obtained through ultrasonic experiments. The longitudinal wave v of ultrasonic waves propagating in the medium P and the shear wave velocity v s , and then the macroscopic bulk modulus K of the sample can be estimated by the following formula: sample and shear modulus μ sample
[0101]
[0102]
[0103] The objective function of constructing numerical simulation data and ultrasonic experimental data based on the least squares method is shown in the following formula:
[0104]
[0105] Where s represents the input parameter data set to be calibrated, including the final rock skeleton V matrix The micro spring stiffness parameter, that is, the normal spring stiffness k n and the tangential spring stiffness k s , K(s) and μ(s) represent the simulation results obtained with different s. After multiple iterations of calibration, the optimal input parameter data set s is obtained, that is, the lattice spring model for numerical simulation of seismic response of organic-rich shale taking into account the microstructure of the organic clay complex is obtained.
[0106] The present invention has been described in detail above. It will be apparent to those skilled in the art that the present invention may be implemented in a wide range under equivalent parameters, concentrations and conditions without departing from the spirit and scope of the present invention and without the need for unnecessary experimentation. Although the present invention provides specific embodiments, it should be understood that further improvements may be made to the present invention. In short, according to the principles of the present invention, this application intends to include any changes, uses or improvements to the present invention, including changes made by conventional techniques known in the art that depart from the scope disclosed in this application. Applications of some of the basic features may be made within the scope of the following appended claims.
Claims
1. A numerical simulation method for seismic response of organic-rich shale, characterized in that: The following steps are involved: 1) Through information collection and image processing of organic-rich shale core samples, digital core two-dimensional images of core samples are obtained and the main mineral components and spatial distribution characteristics of core samples are captured; 2) Generate a rock physics model of the organic-rich shale sample based on the main mineral components and spatial distribution characteristics of the core sample obtained in step 1; 3) Establishing a corresponding lattice spring model for seismic response simulation modeling according to the rock physics model generated in step 2; The step 2) comprises: 21) Statistical analysis of the distribution characteristics of organic matter, organic clay complexes, and clay minerals in core samples, determine a reasonable random medium statistical model, and generate a three-dimensional random medium image based on the size of the digital core sample, denoted as V sample , using a mixed random medium three-dimensional model: Where a, b, and c are the autocorrelation lengths of the medium in the x, y, and z directions, and r is called the roughness factor. These parameters are estimated by inverting the scanned image; 22) Perform Boolean operations on three-dimensional random medium images, V rest =V sample -(V pore ∪V pyrite ) (3) The three-dimensional image after Boolean operation is the preliminary rock skeleton, which is denoted as V rest Binarized image processing and watershed transformation were performed to divide the preliminary rock skeleton into clusters representing organic matter, viscosity minerals, and organic clay complexes, as well as the final rock skeleton, and the corresponding geometric information vector files were obtained, which were recorded as V1, ..., V i ,……,V N , V matrix , where the subscript i represents the i-th cluster, N represents the number of identified clusters, and V matrix For the final rock skeleton; 23) Calculate the irregular pore structure V that makes up the digital core sample pore , Pyrite V pyrite , clusters V1, ..., V generated based on random media i ,……,V N , and finally the rock skeleton V matrix Elastic material parameters, These parameters satisfy the following formula: V sample =V1+...+V i +...+V N +V matrix +V pore +V pyrite (4) The pore structure V pore With pyrite V pyrite The elastic material parameters of the clump are considered to be composed of one component; the clump and the final rock skeleton are described by two independent material elastic parameters; assuming that the components of the clump are mainly composed of clay and organic matter, the proportions of the two in the i-th clump are recorded as and Both meet When hour, Indicates that the mass is all organic matter. hour, Indicates that the lumps are all pure clay. and When , it means that the agglomerate is an organic clay complex, and each agglomerate has and The elastic material parameters of each agglomerate are estimated based on the average value of the Hashin-Shtrikman boundary theory according to the clay and organic matter content of the sample. The volume model and shear modulus of the i-th agglomerate are expressed as K i With μ i express:
2. The method for numerical simulation of seismic response of organic-rich shale according to claim 1, characterized in that: The step 1) comprises: 11) Collecting information and processing images of the core sample includes converting the digital core two-dimensional image of the core sample into a grayscale image using the following color-to-grayscale value formula (1): I Gray =I Red ·0.299+I Green ·0.587+I Blue ·0.114 (1) In the formula, I Gray Represents gray value; I Red ,I Green ,I Blue Represent the red, green, and blue component values respectively; 12) Determine the main minerals and spatial distribution characteristics of the core samples based on the differences in X-ray transmittance of different media; 13) Reconstructing the digital core two-dimensional image of the core sample into a three-dimensional image by combining the grayscale image obtained in step 11) and the main mineral and spatial distribution characteristics of the core sample obtained in step 12), and naming the three-dimensional image digital core.
3. The method for numerical simulation of seismic response of organic-rich shale according to claim 2, characterized in that: In the step 12), the medium of the core sample is divided into three groups of media: irregular pore structure, pyrite and preliminary rock skeleton.
4. The method for numerical simulation of seismic response of organic-rich shale according to claim 2, characterized in that: The information collection of the core samples includes using high-precision CT imaging, comprehensive mineral analysis electron microscope and / or scanning electron microscope.
5. The method for numerical simulation of seismic response of organic-rich shale according to claim 4, characterized in that: The structural materials described in step 22) include organic matter, viscosity minerals, agglomerates of organic clay complexes and the final rock skeleton.
6. The method for numerical simulation of seismic response of organic-rich shale according to claim 1, characterized in that: The step 3) comprises: 31) Generate randomly and closely arranged lattice nodes according to the size of the digital core, determine the lattice nodes representing each group of minerals, and classify and number them; 32) generating a lattice spring network based on a nearest neighbor search algorithm, and determining input parameters of each lattice node in the lattice spring model; 33) Ultrasonic testing is performed on the core sample in step 1), and the micro spring stiffness parameters representing the final rock skeleton are identified by integrating the ultrasonic experimental data. A numerical simulation of the seismic response is performed based on the lattice spring model, thereby obtaining a lattice spring model for numerical simulation of the seismic response of organic-rich shale taking into account the microstructure of the organic clay complex.
7. The method for numerical simulation of seismic response of organic-rich shale according to any one of claims 1 to 6, characterized in that: The numerical simulation method for the seismic response of organic-rich shale is a numerical simulation method for the seismic response of organic-rich shale that can fully consider the microstructure of the organic clay complex.