本发明涉及油气勘探开发技术领域，具体地指一种基于成岩过程的碎屑岩储层孔隙度预测方法。
背景技术：
储层孔隙度是决定油气藏规模和开采价值的重要物性参数。获取储层孔隙度主要依托以下几种地质资料：
1)岩心资料，在实验室中通过一些分析测试手段对单个岩心进行测量，从而获取岩石的孔隙度参数；
2)测井资料，利用岩石的电性特征来间接获取岩石的孔隙度；
3)地震资料，利用砂岩波阻抗进行砂岩孔隙度的计算。
通过岩心资料和测井资料获取的孔隙度数据常常有限，井点与井点之间的孔隙度数据需进行数学插值，而地下储层非均质性常制约数据的可靠性，另外受采样点的限制，所采集数据不一定都有代表性。地震数据具有很好连续性，而限于品质和分辨率原因，不能满足储层精细研究和刻画，尤其是砂体薄的碎屑岩储层。基于这一背景，利用直接的储层测量数据已不能较好地反映出储层的发育，特别是有效储层的发育分布状况，影响了碎屑岩储层，尤其是高成岩阶段低孔、低渗储层油气勘探开发，必须寻找出一种不同于常规方法的评价方法，来更好地、更合理地预测优势储层分布，最终为油气勘探开发提供有力依据。
成岩作用是一切储层形成必经过程，有什么样的作用过程，就有什么样的响应特征。储层中物性变化是地质历史中压实、胶结以及溶蚀作用等的共同结果。考虑储层演化过程中不同时期储层发生的各种胶结作用，特别是不同时期胶结作用发生时的条件及其差异，以及这些条件及差异作用下的胶结结果，可以为处于碎屑岩储层孔隙预测提供很好的手段。目前国内外储层成岩数值模拟原理与技术存在以下几点不足：
1)在对孔隙进行数值模拟时，地质参数选取单一，没能综合考虑沉积相、粒度大小、岩屑颗粒、刚性颗粒、磨圆度、分选性、泥质含量、胶结强度、胶结类型、胶结时间、溶蚀时间、埋藏温度、异常高压、沉积速率等多种地质参数对孔隙影响；
2)模拟过程过于线性，不能再现储层演化的复杂过程；
3)模拟方法过于综合，不能满足科研需要以及油气藏精细描述。如果不能建立合理地质模型，综合考虑各种地质参数条件下压实作用、胶结作用、溶蚀作用之间影响关系及其主次之分，一起数值模拟就只会是数值游戏，也就不能真正解决实际问题。
技术实现要素：
本发明针对上述存在的问题，提供了基于成岩过程的碎屑岩储层孔隙度预测方法，该方法基于碎屑岩储层沉积相、岩性、流体性质、成岩阶段的综合研究，预测不同成岩阶段下储层与成岩作用时空匹配关系，建立不同地质环境条件下以及不同成岩作用下“地质参数-成岩作用-孔隙度”预测模型，进而对碎屑岩储层孔隙度大小进行预测，确定碎屑岩储层孔隙空间分布，为油藏评价提供依据
为了实现上述目的，本发明提供了一种基于成岩过程的碎屑岩储层孔隙度预测方法，包括以下步骤：
1)收集研究区地质资料
2)建立原始资料数据库，包括
(1)沉积相类型数据库f_m
对碎屑岩储层沉积相类型进行概括和分类建立沉积相类型数据库f_m
(2)岩性数据库r_n
根据组成岩石颗粒的颗粒大小及矿物成分，将碎屑岩岩性进行划分，建立岩性数据库，
(3)流体数据库p_o
根据储层孔隙中流体性质，建立流体数据库；
(4)成岩阶段数据库s_q
(5)建立成岩作用数据库d_c
根据岩石中不同矿物mi的孔隙改变量建立成岩作用数据库；
(6)建立成岩作用数学模型
成岩作用数学模型包括
a.压实级别数学模型a_k
b.胶结级别数学模型b_k
c.溶蚀级别数学模型c_k
3)建立成岩阶段预测模型
(1)确定演化期次
根据目的层上覆地层发育情况划分目标层在地质历史时期所经历的演化阶段，目的层l上覆地层有i层，即从上到下依次标记为：l1、l2、l3……li-1、li，则目的层在地质历史时期演化阶段共有i个，按照时间演化，目的层l演化期次先后顺序依次计为li、li-1、li-2……l2、l1；
(2)确定不同期次埋深
在(1)的基础下，计算目的层l不同演化期次的埋深，其计算公式如下：
dep(li)＝h(l)-h(li)；
dep(li-1)＝h(l)-h(li-1)；
dep(li-2)＝h(l)-h(li-2)；
……
dep(l2)＝h(l)-h2；
dep(l1)＝h(l)-h1；
其中：h(l)为目的层l顶界面，h(li)为上覆层li顶界面，dep(li)为目的层li阶段埋深；
(3)确定不同期次地层温度
地质体在埋藏的过程中，温度的大小可表示为与深度的线性函数关系，通过该模型，可计算目的层l在不同时期、不同深度、不同位置地层温度；其温度计算模型公式：
t＝t0+c*(d(ti)-h0)
其中t0为常温带温度，c为常数，dep(ti)为目的层ti阶段埋深，h0为恒温带埋深，为常数，t为目标层古地温；
(4)确定不同期次地层成岩阶段
4)建立成岩作用预测模型
(1)研究区储层网格化
将储层网格化，即空间上不均匀分布的数据按一定方法(如滑动平均法、克里格法或其他适当的数值推算方法)归算成规则网格中的代表值(趋势值)的过程；研究区储层的每个网格可用wi(x，y)表示；
(2)确定网格属性
a.根据步骤3)确定网格wi(x，y)不同演化期次成岩阶段s_q；
b.根据研究区沉积相数据确定网格wi(x，y)的沉积相属性f_m；
c.根据研究区岩性数据确定网格wi(x，y)的岩性属性r_n；
d.根据研究区流体数据确定网格wi(x，y)的流体属性p_o；
e.根据上述步骤a、b、c和d确定每个网格不同演化期次wi(x，y)的综合属性g(s_q，f_m，r_n，p_o)，即：wi(x，y)＝s_q+f_m+r_n+p_o；
(3)基于研究区地质资料，确定研究区已知井hj所在网格wj(x，y)不同演化期次的成岩阶段s_q_j以及成岩作用d_s_j；
(4)确定已知井hj所在网格wj(x，y)不同演化期次不同演化阶段的综合属性gj(f_m_j，r_n_j，p_o_j)，建立已知井hj所在网格wj(x，y)不同演化期次不同演化阶段成岩作用d_s与综合属性g(f_m_j，r_n_j，p_o_j)的对应关系，即为：
wj(s_q_j，d_s_j)＝f_m_j+r_n_j+p_o_j；
(5)任取一未知网格wi(x，y)，确定未知网格wi(x，y)成岩阶段s_q时的综合属性g(f_m，r_n，p_o)，即为：
gi＝(f_m_i，r_n_i，p_o_i)；
(6)将网格wi(x，y)成岩阶段s_q时的成岩作用综合属性gi与井hj所在网格wj(x，y)成岩阶段s_q的综合属性gj相比较，即为：
gi-gj＝(f_m_i，r_n_i，p_o_i)-(f_m_j+r_n_j+p_o_j)；
如果满足gi-gj＝0，则未知网格wi(x，y)的成岩作用d_s_i与井j具有相同的成岩作用d_s_j；
如果gi-gj≠0，则按不同属性优先级顺序，即为：沉积相f_m一级、岩性r_n二级、流体性质p_o三级对未知网格wi(x，y)的成岩作用进行判识，即：
ⅰ：f_m_i-f_m_j＝0，r_n_i-r_n_j≠0，p_o_i-p_o_j≠0；
ⅱ：f_m_i-f_m_j＝0，r_n_i-r_n_j＝0，p_o_i-p_o_j≠0；
未知网格wi(x，y)的成岩作用d_s_i满足条件i时，未知网格wi(x，y)与具有相同属性f_j的网格wj(x，y)的成岩作用相同；
d_s_i满足条件ⅱ时，未知网格wi(x，y)与具有相同属性f_m_j、r_n_j的网格wj(x，y)的成岩作用相同；
(7)根据每个网格wi(x，y)不同演化期次的成岩作用确定研究区成岩作用演化；
5)计算不同演化期次压实量、胶结量和溶蚀量
(1)根据成岩序列数据，确定已知井所在网格wj(x，y)的成岩阶段s_q所对应的压实作用、胶结作用、溶蚀作用下的压实级别数学模型a_k_j、胶结级别数学模型b_k_j、溶蚀级别数学模型c_k_j，即wj(s_q_j，d_s_j)＝wj(a_k_j，b_k_j，c_k_j)；
(2)将wi(x，y)的wi(s_q_j，d_s_j)与wj(s_q_j，d_s_j)比对，满足wi(s_q_j，d_s_j)＝wj(s_q_j，d_s_j)，则wi(x，y)与wj(s_q_j，d_s_j)具有相同的压实级别数学模型a_k、胶结级别数学模型b_k、溶蚀级别数学模型c_k，即a_k_i＝a_k_j、b_k_i＝b_k_j、c_k_i＝c_k_j；如果出现同一未知网格wi(x，y)与多个已知网格wj(s_q_j，d_s_j)相同，则a_k_i＝1/n∑a_k_j、b_k_i＝1/n∑b_k_j、c_k_i＝1/n∑c_k_j，n为相同网格wj(s_q_j，d_s_j)的网格数；
(3)根据网格wi(x，y)不同演化期次的压实级别数学模型a_k、胶结级别数学模型b_k、溶蚀级别数学模型c_k计算每个网格wi(x，y)不同演化期次的压实率rcom、胶结率rmi_cem、溶蚀率rmi_dis(mi表示石英、长石、方解石、白云石、黏土矿物)；
(4)根据每个网格wi(x，y)的压实率rcom、胶结率rmi_cem和溶蚀率rmi_dis计算压实量vcom，胶结量vmi_cem，溶蚀量vmi_dis，
其中，
(5)根据每个网格wi(x，y)不同演化期次的孔隙值确定目标储层孔隙值空间分布。
进一步地，所述步骤1)中，收集研究区地质资料包括：
研究区构造等值线图、研究区沉积相图、研究区岩性分布图、研究区流体分布图。
再进一步地，所述步骤2)中，
(1)沉积相类型数据库f_m包括冲积扇-旱扇-扇根f_1、冲积扇-旱扇-扇中f_2、冲积扇-旱扇-扇缘f_3、冲积扇-湿扇-扇根f_4、冲积扇-湿扇-扇中f_5、冲积扇-湿扇-扇缘f_6、河流相-曲流河-河床亚相-河床滞留f_7、河流相-曲流河-河床亚相-边滩f_8、河流相-曲流河-堤岸亚相-天然堤f_9、河流相-曲流河-堤岸亚相-决口扇f_10、河流相-曲流河-河漫亚相-河漫滩f_11、河流相-曲流河-泛滥盆地-河漫湖泊f_12、河流相-曲流河-河漫沼泽f_13、河流相-辫状河-牛轭湖f_14、河流相-辫状河-河床滞留f_15、河流相-辫状河-心滩f_16、河道f_17、河流相-辫状河-泛滥平原f_18、湖泊相-断陷型-湖成三角洲f_19、湖泊相-断陷型-滨湖f_20、湖泊相-断陷型-浅湖f_21、湖泊相-断陷型-半深湖f_22、湖泊相-断陷型-深湖f_23、湖泊相-断陷型-湖湾f_24、湖泊相-坳陷型-湖成三角洲f_25、湖泊相-坳陷型-滨湖f_26、湖泊相-坳陷型-浅湖f_27、湖泊相-坳陷型-半深湖f_28、湖泊相-坳陷型-深湖f_29、湖泊相-坳陷型-湖湾f_30、湖泊相-前陆型-湖成三角洲f_31、湖泊相-前陆型-滨湖f_32、湖泊相-前陆型-浅湖f_33、湖泊相-前陆型-半深湖f_34、湖泊相-前陆型-深湖f_35、湖泊相-前陆型-湖湾f_36、三角洲相-辫状河三角洲-三角洲平原-分支f_37、三角洲相-辫状河三角洲-三角洲平原-河道f_38、三角洲相-辫状河三角洲-三角洲平原-陆上天然堤f_39、三角洲相-辫状河三角洲-三角洲平原-决口扇f_40、三角洲相-辫状河三角洲-三角洲平原-沼泽f_41、三角洲相-辫状河三角洲-三角洲平原-淡水湖泊f_42、三角洲相-辫状河三角洲-三角洲前缘-水下分支河道f_43、三角洲相-辫状河三角洲-三角洲前缘-水下天然堤f_44、三角洲相-辫状河三角洲-三角洲前缘-支流间湾f_45、三角洲相-辫状河三角洲-三角洲前缘-分支河口砂坝f_46、三角洲相-辫状河三角洲-三角洲前缘-远砂坝f_47、三角洲相-辫状河三角洲-三角洲前缘-前缘席状砂f_48、三角洲相-辫状河三角洲-前三角洲f_49、三角洲相-扇三角洲-扇三角洲平原-泥石流f_50、三角洲相-扇三角洲-扇三角洲平原-河道充填f_51、三角洲相-扇三角洲-扇三角洲平原-漫滩f_52、三角洲相-扇三角洲-扇三角洲前缘-碎屑流f_53、三角洲相-扇三角洲-扇三角洲前缘-水下分流河道f_54、三角洲相-扇三角洲-扇三角洲前缘-支流间湾f_55、三角洲相-扇三角洲-扇三角洲前缘-河口砂坝f_56、三角洲相-扇三角洲-扇三角洲前缘-远砂坝f_57、三角洲相-扇三角洲-前扇三角洲f_58；
(2)岩性数据库r_n包括细砾岩r_1、中砾岩r_2、粗砾岩r_3、巨砾岩r_4、石英砂岩r_5、长石质石英砂岩r_6、岩屑质石英砂岩r_7、长石岩屑质石英砂岩r_8、长石砂岩r_9、岩屑质长石砂岩r_10、岩屑长石砂岩r_11、岩屑砂岩r_12、长石质岩屑砂岩r_13、长石岩屑砂岩r_14、黏土r_15、泥岩r_16、页岩r_17；
(3)流体数据库p_o包括淡水环境p_1、酸性水环境p_2、碱性水环境p_3；
(4)成岩阶段数据库s_q包括：早成岩阶段a期ⅰa，即为：s_1；早成岩阶段b期ⅰb，即为：s_2；中成岩阶段a期ⅱa，即为：s_3；中成岩阶段b期ⅱb，即为s_4；晚成岩阶段ⅲ，即为：s_5；
(5)建立成岩作用数据库d_c
根据岩石中不同矿物mi的孔隙改变量建立成岩作用数据库，成岩作用类型包括：
强压实作用d1-1：其孔隙改变量为
中压实作用d1-2：其孔隙改变量为
弱压实作用d1-3：其孔隙改变量
石英强胶结作用d2-1：其孔隙改变量为
石英中胶结作用d2-2：其孔隙改变量为
石英弱胶结作用d2-3：其孔隙改变量
长石强胶结作用d3-1：孔隙改变量
长石中胶结作用d3-2：其孔隙改变量
长石弱胶结作用d3-3：其孔隙改变量
方解石强胶结作用d4-1：其孔隙改变量
方解石中胶结作用d4-2：其孔隙改变量
方解石弱胶结作用d4-3：其孔隙改变量
白云石强胶结作用d5-1：其孔隙改变量
白云石中胶结作用d5-2：其孔隙改变量
白云石弱胶结作用d5-3：其孔隙改变量
黏土强胶结作用d6-1：其孔隙改变量
黏土中胶结作用d6-2：其孔隙改变量
黏土弱胶结作用d6-3：其孔隙改变量
石英强溶蚀作用d7-1：其孔隙改变量为
石英中溶蚀作用d7-2：其孔隙改变量为
石英弱溶蚀作用d7-3：其孔隙改变量
长石强溶蚀作用d8-1：孔隙改变量
长石中溶蚀作用d8-2：其孔隙改变量
长石弱溶蚀作用d8-3：其孔隙改变量
方解石强溶蚀作用d9-1：其孔隙改变量
方解石中溶蚀作用d9-2：其孔隙改变量
方解石弱溶蚀作用d9-3：其孔隙改变量
白云石强溶蚀作用d10-1：其孔隙改变量
白云石中溶蚀作用d10-2：其孔隙改变量
白云石弱溶蚀作用d10-3：其孔隙改变量
(6)建立成岩作用数学模型
a.压实级别数学模型a_k
a_1：ⅰ：y＝1.9467*ln(dep)-0.8689；
a_2：ⅱ：y＝2.9200*ln(dep)-1.3034；
a_3：ⅲ：y＝4.3800*ln(dep)-1.9551；
a_4：ⅳ：y＝7.0080*ln(dep)-3.1282；
a_5：ⅴ：y＝8.7600*ln(dep)-3.9103；
a_6：ⅵ：y＝11.6800*ln(dep)-5.2137；
a_7：ⅶ：y＝14.0160*ln(dep)-6.2564；
a_8：ⅷ：y＝17.5200*ln(dep)-7.8205；
a_9：ⅸ：y＝20.6118*ln(dep)-9.2006；
a_10：ⅹ：y＝23.3600*ln(dep)-10.4273；
其中，y压实率，ma，y胶结率，％；
b.胶结级别数学模型b_k
建立起不同类型胶结物mi的数学模型，从而定量刻画碎屑岩在成岩过程中孔隙减小量；
b_1＝ⅰ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/10；
b_2＝ⅱ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/8；
b_3＝ⅲ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/6；
b_4＝ⅳ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/5；
b_5＝ⅴ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/3；
b_6＝ⅵ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.6；
b_7＝ⅶ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.5；
b_8＝ⅷ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.4；
b_9＝ⅸ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.2；
b_10＝ⅹ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1；
其中，x埋藏时间，ma，y胶结率，％；
c.溶蚀级别数学模型c_k
建立起不同类型溶蚀物mi的数学模型，从而定量刻画碎屑岩在成岩过程中孔隙减小量；
c_1＝ⅰ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/10；
c_2＝ⅱ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/8；
c_3＝ⅲ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/6；
c_4＝ⅳ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/5；
c_5＝ⅴ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/3；
c_6＝ⅵ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.6；
c_7＝ⅶ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.5；
c_8＝ⅷ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.4；
c_9＝ⅸ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.2；
c_10＝ⅹ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1；
其中，x埋藏时间，ma，y溶蚀率，％。
再进一步地，所述步骤2)第(4)小步中，不同期次地层成岩阶段根据不同温度确定：
当t∈[20～65)，ds为早成岩阶段早期，即为ⅰa；
当t∈[65～85)，ds为早成岩阶段晚期，即为ⅰb；
当t∈[85～140)，ds为中成岩阶段早期，即ⅱa；
当t∈[140～175)，ds为中成岩阶段晚期，即为ⅱb；
当t∈[175～200)，ds为晚成岩阶段，即为ⅲ。
本发明的有益效果在于：
本发明提供了基于成岩过程的碎屑岩储层孔隙度预测方法，该方法基于碎屑岩储层沉积相、岩性、流体性质、成岩阶段的综合研究，预测不同成岩阶段下储层与成岩作用时空匹配关系，建立不同地质环境条件下以及不同成岩作用下“地质参数-成岩作用-孔隙度”预测模型，进而对碎屑岩储层孔隙度大小进行预测，确定碎屑岩储层孔隙空间分布。该方法为一种不同于常规方法的评价方法，能更好地、更合理地预测优势储层分布，最终为油气勘探开发提供有力依据。
附图说明
图1为研究区平原组构造等值线图；
图2为研究区明化镇组构造等值线图；
图3为研究区馆陶组构造等值线图；
图4为研究区东营组二段构造等值线图；
图5为研究区东营组二段沉积相图；
图6为研究区东营组二段岩性分布图；
图7为研究区馆陶组沉积时东营组二段流体分布图；
图8为研究区明化镇组沉积时东营组二段流体分布图；
图9为研究区平原组沉积时东营组二段流体分布图；
图10为研究区馆陶组沉积时东营组二段地温图；
图11为研究区明化镇组沉积时东营组二段地温图；
图12为研究区平原组沉积时东营组二段地温图；
图13为研究区馆陶组沉积时东营组二段成岩阶段图；
图14为研究区明化镇组沉积时东营组二段成岩阶段图；
图15为研究区平原组沉积时东营组二段成岩阶段图；
图16为研究区馆陶组沉积时东营组二段成岩作用图；
图17为研究区明化镇组沉积时东营组二段成岩作用图；
图18为研究区平原组沉积时东营组二段成岩作用图；
图19为研究区馆陶组沉积时东营组二段孔隙度分布图；
图20为研究区明化镇组沉积时东营组二段孔隙度分布图；
图21为研究区平原组沉积时东营组二段孔隙度分布图。
具体实施方式
为了更好地解释本发明，以下结合具体实施例进一步阐明本发明的主要内容，但本发明的内容不仅仅局限于以下实施例。
实施例1
基于成岩过程的碎屑岩储层孔隙度预测方法，包括以下步骤：
1)收集研究区地质资料，它包括：
①研究区东营组二段压实作用表、胶结作用表、溶蚀作用表；
②研究区平原组构造等值线图；
③研究区明化镇组构造等值线图；
④研究区馆陶组构造等值线图；
⑤研究区东营组二段构造等值线图；
⑥研究区东营组二段沉积相图；
⑦研究区东营组二段岩性分布图；
⑧研究区馆陶组沉积时东营组二段流体分布图；
⑨研究区明化镇组沉积时东营组二段流体分布图；
⑩研究区平原组沉积时东营组二段流体分布图；
2)建立原始资料数据库
(1)沉积相类型数据库f_m
对碎屑岩储层沉积相类型进行概括和分类，其中包括：冲积扇-旱扇-扇根f_1、冲积扇-旱扇-扇中f_2、冲积扇-旱扇-扇缘f_3、冲积扇-湿扇-扇根f_4、冲积扇-湿扇-扇中f_5、冲积扇-湿扇-扇缘f_6、河流相-曲流河-河床亚相-河床滞留f_7、河流相-曲流河-河床亚相-边滩f_8、河流相-曲流河-堤岸亚相-天然堤f_9、河流相-曲流河-堤岸亚相-决口扇f_10、河流相-曲流河-河漫亚相-河漫滩f_11、河流相-曲流河-泛滥盆地-河漫湖泊f_12、河流相-曲流河-河漫沼泽f_13、河流相-辫状河-牛轭湖f_14、河流相-辫状河-河床滞留f_15、河流相-辫状河-心滩f_16、河道f_17、河流相-辫状河-泛滥平原f_18、湖泊相-断陷型-湖成三角洲f_19、湖泊相-断陷型-滨湖f_20、湖泊相-断陷型-浅湖f_21、湖泊相-断陷型-半深湖f_22、湖泊相-断陷型-深湖f_23、湖泊相-断陷型-湖湾f_24、湖泊相-坳陷型-湖成三角洲f_25、湖泊相-坳陷型-滨湖f_26、湖泊相-坳陷型-浅湖f_27、湖泊相-坳陷型-半深湖f_28、湖泊相-坳陷型-深湖f_29、湖泊相-坳陷型-湖湾f_30、湖泊相-前陆型-湖成三角洲f_31、湖泊相-前陆型-滨湖f_32、湖泊相-前陆型-浅湖f_33、湖泊相-前陆型-半深湖f_34、湖泊相-前陆型-深湖f_35、湖泊相-前陆型-湖湾f_36、三角洲相-辫状河三角洲-三角洲平原-分支f_37、三角洲相-辫状河三角洲-三角洲平原-河道f_38、三角洲相-辫状河三角洲-三角洲平原-陆上天然堤f_39、三角洲相-辫状河三角洲-三角洲平原-决口扇f_40、三角洲相-辫状河三角洲-三角洲平原-沼泽f_41、三角洲相-辫状河三角洲-三角洲平原-淡水湖泊f_42、三角洲相-辫状河三角洲-三角洲前缘-水下分支河道f_43、三角洲相-辫状河三角洲-三角洲前缘-水下天然堤f_44、三角洲相-辫状河三角洲-三角洲前缘-支流间湾f_45、三角洲相-辫状河三角洲-三角洲前缘-分支河口砂坝f_46、三角洲相-辫状河三角洲-三角洲前缘-远砂坝f_47、三角洲相-辫状河三角洲-三角洲前缘-前缘席状砂f_48、三角洲相-辫状河三角洲-前三角洲f_49、三角洲相-扇三角洲-扇三角洲平原-泥石流f_50、三角洲相-扇三角洲-扇三角洲平原-河道充填f_51、三角洲相-扇三角洲-扇三角洲平原-漫滩f_52、三角洲相-扇三角洲-扇三角洲前缘-碎屑流f_53、三角洲相-扇三角洲-扇三角洲前缘-水下分流河道f_54、三角洲相-扇三角洲-扇三角洲前缘-支流间湾f_55、三角洲相-扇三角洲-扇三角洲前缘-河口砂坝f_56、三角洲相-扇三角洲-扇三角洲前缘-远砂坝f_57、三角洲相-扇三角洲-前扇三角洲f_58；
(2)岩性数据库r_n
根据组成岩石颗粒的颗粒大小及矿物成分，将碎屑岩岩性进行划分，建立岩性数据库，它包括细砾岩r_1、中砾岩r_2、粗砾岩r_3、巨砾岩r_4、石英砂岩r_5、长石质石英砂岩r_6、岩屑质石英砂岩r_7、长石岩屑质石英砂岩r_8、长石砂岩r_9、岩屑质长石砂岩r_10、岩屑长石砂岩r_11、岩屑砂岩r_12、长石质岩屑砂岩r_13、长石岩屑砂岩r_14、黏土r_15、泥岩r_16、页岩r_17；
(3)流体数据库p_o
根据储层孔隙中流体性质，建立流体数据库：它包括淡水环境p_1、酸性水环境p_2、碱性水环境p_3；
(4)成岩阶段数据库s_q
成岩阶段数据库s_q包括：早成岩阶段a期ⅰa，即为：s_1；早成岩阶段b期ⅰb，即为：s_2；中成岩阶段a期ⅱa，即为：s_3；中成岩阶段b期ⅱb，即为s_4；晚成岩阶段ⅲ，即为：s_5；
(5)建立成岩作用数据库d_c
根据岩石中不同矿物mi的孔隙改变量建立成岩作用数据库，成岩作用类型包括：
强压实作用d1-1：其孔隙改变量为
中压实作用d1-2：其孔隙改变量为
弱压实作用d1-3：其孔隙改变量
石英强胶结作用d2-1：其孔隙改变量为
石英中胶结作用d2-2：其孔隙改变量为
石英弱胶结作用d2-3：其孔隙改变量
长石强胶结作用d3-1：孔隙改变量
长石中胶结作用d3-2：其孔隙改变量
长石弱胶结作用d3-3：其孔隙改变量
方解石强胶结作用d4-1：其孔隙改变量
方解石中胶结作用d4-2：其孔隙改变量
方解石弱胶结作用d4-3：其孔隙改变量
白云石强胶结作用d5-1：其孔隙改变量
白云石中胶结作用d5-2：其孔隙改变量
白云石弱胶结作用d5-3：其孔隙改变量
黏土强胶结作用d6-1：其孔隙改变量
黏土中胶结作用d6-2：其孔隙改变量
黏土弱胶结作用d6-3：其孔隙改变量
石英强溶蚀作用d7-1：其孔隙改变量为
石英中溶蚀作用d7-2：其孔隙改变量为
石英弱溶蚀作用d7-3：其孔隙改变量
长石强溶蚀作用d8-1：孔隙改变量
长石中溶蚀作用d8-2：其孔隙改变量
长石弱溶蚀作用d8-3：其孔隙改变量
方解石强溶蚀作用d9-1：其孔隙改变量
方解石中溶蚀作用d9-2：其孔隙改变量
方解石弱溶蚀作用d9-3：其孔隙改变量
白云石强溶蚀作用d10-1：其孔隙改变量
白云石中溶蚀作用d10-2：其孔隙改变量
白云石弱溶蚀作用d10-3：其孔隙改变量
(6)建立成岩作用数学模型
a.压实级别数学模型a_k
a_1：ⅰ：y＝1.9467*ln(dep)-0.8689；
a_2：ⅱ：y＝2.9200*ln(dep)-1.3034；
a_3：ⅲ：y＝4.3800*ln(dep)-1.9551；
a_4：ⅳ：y＝7.0080*ln(dep)-3.1282；
a_5：ⅴ：y＝8.7600*ln(dep)-3.9103；
a_6：ⅵ：y＝11.6800*ln(dep)-5.2137；
a_7：ⅶ：y＝14.0160*ln(dep)-6.2564；
a_8：ⅷ：y＝17.5200*ln(dep)-7.8205；
a_9：ⅸ：y＝20.6118*ln(dep)-9.2006；
a_10：ⅹ：y＝23.3600*ln(dep)-10.4273；
其中，y压实率，ma，y胶结率，％；
b.胶结级别数学模型b_k
建立起不同类型胶结物mi的数学模型，从而定量刻画碎屑岩在成岩过程中孔隙减小量；
b_1＝ⅰ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/10；
b_2＝ⅱ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/8；
b_3＝ⅲ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/6；
b_4＝ⅳ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/5；
b_5＝ⅴ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/3；
b_6＝ⅵ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.6；
b_7＝ⅶ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.5；
b_8＝ⅷ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.4；
b_9＝ⅸ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1.2；
b_10＝ⅹ：y＝[7.515e-05*(x)^3-0.0158*(x)^2+1.686*(x)+14.38]/1；
其中，x埋藏时间，ma，y胶结率，％；
c.溶蚀级别数学模型c_k
建立起不同类型溶蚀物mi的数学模型，从而定量刻画碎屑岩在成岩过程中孔隙减小量；
c_1＝ⅰ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/10；
c_2＝ⅱ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/8；
c_3＝ⅲ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/6；
c_4＝ⅳ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/5；
c_5＝ⅴ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/3；
c_6＝ⅵ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.6；
c_7＝ⅶ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.5；
c_8＝ⅷ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.4；
c_9＝ⅸ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1.2；
c_10＝ⅹ：y＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/1；
其中，x埋藏时间，ma，y溶蚀率，％；
3)建立成岩阶段预测模型
(1)确定演化期次
根据目的层上覆地层发育情况划分目标层在地质历史时期所经历的演化阶段，目的层l上覆地层有i层，即从上到下依次标记为：l1、l2、l3……li-1、li，则目的层在地质历史时期演化阶段共有i个，按照时间演化，目的层l演化期次先后顺序依次计为li、li-1、li-2……l2、l1；
(2)确定不同期次埋深
在(1)的基础下，计算目的层l不同演化期次的埋深，其计算公式如下：
dep(li)＝h(l)-h(li)；
dep(li-1)＝h(l)-h(li-1)；
dep(li-2)＝h(l)-h(li-2)；
……
dep(l2)＝h(l)-h2；
dep(l1)＝h(l)-h1；
其中：h(l)为目的层l顶界面，h(li)为上覆层li顶界面，dep(li)为目的层li阶段埋深；
(3)确定不同期次地层温度
地质体在埋藏的过程中，温度的大小可表示为与深度的线性函数关系，通过该模型，可计算目的层l在不同时期、不同深度、不同位置地层温度；其温度计算模型公式：
t＝t0+c*(d(ti)-h0)
其中t0为常温带温度，c为常数，dep(ti)为目的层ti阶段埋深，h0为恒温带埋深，为常数，t为目标层古地温；
(4)确定不同期次地层成岩阶段
当t∈[20～65)，ds为早成岩阶段早期，即为ⅰa；
当t∈[65～85)，ds为早成岩阶段晚期，即为ⅰb；
当t∈[85～140)，ds为中成岩阶段早期，即ⅱa；
当t∈[140～175)，ds为中成岩阶段晚期，即为ⅱb；
当t∈[175～200)，ds为晚成岩阶段，即为ⅲ；
4)建立成岩作用预测模型
(1)研究区储层网格化
将储层网格化，即空间上不均匀分布的数据按一定方法(如滑动平均法、克里格法或其他适当的数值推算方法)归算成规则网格中的代表值(趋势值)的过程；研究区储层的每个网格可用wi(x，y)表示；(2)确定网格属性
a.根据步骤3)确定网格wi(x，y)不同演化期次成岩阶段s_q；
b.根据研究区沉积相数据确定网格wi(x，y)的沉积相属性f_m；
c.根据研究区岩性数据确定网格wi(x，y)的岩性属性r_n；
d.根据研究区流体数据确定网格wi(x，y)的流体属性p_o；
e.根据上述步骤a、b、c和d确定每个网格不同演化期次wi(x，y)的综合属性g(s_q，f_m，r_n，p_o)，即：wi(x，y)＝s_q+f_m+r_n+p_o；
(3)基于研究区地质资料，确定研究区已知井hj所在网格wj(x，y)不同演化期次的成岩阶段s_q_j以及成岩作用d_s_j；
(4)确定已知井hj所在网格wj(x，y)不同演化期次不同演化阶段的综合属性gj(f_m_j，r_n_j，p_o_j)，建立已知井hj所在网格wj(x，y)不同演化期次不同演化阶段成岩作用d_s与综合属性g(f_m_j，r_n_j，p_o_j)的对应关系，即为：
wj(s_q_j，d_s_j)＝f_m_j+r_n_j+p_o_j；
(5)任取一未知网格wi(x，y)，确定未知网格wi(x，y)成岩阶段s_q时的综合属性g(f_m，r_n，p_o)，即为：
gi＝(f_m_i，r_n_i，p_o_i)；
(6)将网格wi(x，y)成岩阶段s_q时的成岩作用综合属性gi与井hj所在网格wj(x，y)成岩阶段s_q的综合属性gj相比较，即为：
gi-gj＝(f_m_i，r_n_i，p_o_i)-(f_m_j+r_n_j+p_o_j)；
如果满足gi-gj＝0，则未知网格wi(x，y)的成岩作用d_s_i与井j具有相同的成岩作用d_s_j；
如果gi-gj≠0，则按不同属性优先级顺序，即为：沉积相f_m一级、岩性r_n二级、流体性质p_o三级对未知网格wi(x，y)的成岩作用进行判识，即：
ⅰ：f_m_i-f_m_j＝0，r_n_i-r_n_j≠0，p_o_i-p_o_j≠0；
ⅱ：f_m_i-f_m_j＝0，r_n_i-r_n_j＝0，p_o_i-p_o_j≠0；
未知网格wi(x，y)的成岩作用d_s_i满足条件i时，未知网格wi(x，y)与具有相同属性f_j的网格wj(x，y)的成岩作用相同；
d_s_i满足条件ⅱ时，未知网格wi(x，y)与具有相同属性f_m_j、r_n_j的网格wj(x，y)的成岩作用相同；
(7)根据每个网格wi(x，y)不同演化期次的成岩作用确定研究区成岩作用演化；
5)计算不同演化期次压实量、胶结量和溶蚀量
(1)根据成岩序列数据，确定已知井所在网格wj(x，y)的成岩阶段s_q所对应的压实作用、胶结作用、溶蚀作用下的压实级别数学模型a_k_j、胶结级别数学模型b_k_j、溶蚀级别数学模型c_k_j，即wj(s_q_j，d_s_j)＝wj(a_k_j，b_k_j，c_k_j)；
(2)将wi(x，y)的wi(s_q_j，d_s_j)与wj(s_q_j，d_s_j)比对，满足wi(s_q_j，d_s_j)＝wj(s_q_j，d_s_j)，则wi(x，y)与wj(s_q_j，d_s_j)具有相同的压实级别数学模型a_k、胶结级别数学模型b_k、溶蚀级别数学模型c_k，即a_k_i＝a_k_j、b_k_i＝b_k_j、c_k_i＝c_k_j；如果出现同一未知网格wi(x，y)与多个已知网格wj(s_q_j，d_s_j)相同，则a_k_i＝1/n∑a_k_j、b_k_i＝1/n∑b_k_j、c_k_i＝1/n∑c_k_j，n为相同网格wj(s_q_j，d_s_j)的网格数；
(3)根据网格wi(x，y)不同演化期次的压实级别数学模型a_k、胶结级别数学模型b_k、溶蚀级别数学模型c_k计算每个网格wi(x，y)不同演化期次的压实率rcom、胶结率rmi_cem、溶蚀率rmi_dis(mi表示石英、长石、方解石、白云石、黏土矿物)；
(4)根据每个网格wi(x，y)的压实率rcom、胶结率rmi_cem和溶蚀率rmi_dis计算压实量vcom，胶结量vmi_cem，溶蚀量vmi_dis，
其中，
(5)根据每个网格wi(x，y)不同演化期次的孔隙值确定目标储层孔隙值空间分布。
实施例2
基于成岩过程的碎屑岩储层孔隙度预测方法预测渤中凹陷东营组二段储层，具体步骤如下：
1)收集研究区地质资料
收集渤中凹陷东营组二段储层的地质资料，并按照原始资料数据库对各种类型数据进行分类。收集的地质资料包括以下内容：研究区东营组二段溶蚀作用表(表1、表2、表3)、研究区平原组构造等值线图(图1)、研究区明化镇组构造等值线图(图2)、研究区馆陶组构造等值线图(图3)、研究区东营组二段构造等值线图(图4)、研究区东营组二段沉积相图(图5)、研究区东营组二段岩性分布图(图6)、研究区馆陶组沉积时东营组二段流体分布图(图7)、研究区明化镇组沉积时东营组二段流体分布图(图8)、研究区平原组沉积时东营组二段流体分布图(图9)。
2)确定成岩阶段
(1)确定演化期次
东营组二段上覆地层从上到下依次是平原组、明化镇组以及馆陶组，即从上到下依次标记为l1、l2、l3，则研究区演化期次共有3个，按时间演化，东营组二段l的演化期次先后顺序依次计为l3、l2、l1；
(2)确定不同期次埋深
在(1)的基础下，计算东营组二段l不同演化期次的埋深，其计算公式如下：
dep(li)＝h(l)-h(li)；
dep(li-1)＝h(l)-h(li-1)；
dep(li-2)＝h(l)-h(li-2)；
……
dep(l2)＝h(l)-h2；
dep(l1)＝h(l)-h1；
其中：h(l)为目的层l顶界面，h(li)为上覆层li顶界面，dep(li)为目的层li阶段埋深；
(3)确定不同期次地层温度
地质体在埋藏的过程中，温度的大小可表示为与深度的线性函数关系，通过该模型，可计算目的层l在不同时期、不同深度、不同位置地层温度；
温度计算模型公式：
t＝t0+c*(d(ti)-h0)
其中t0＝21℃，c＝0.0031℃/m，dep(ti)为目的层ti阶段埋深，h0＝30m，t为目标层古地温；
(4)确定不同期次地层成岩阶段
当t∈[20～65)，ds为早成岩阶段早期，即为ⅰa；
当t∈[65～85)，ds为早成岩阶段晚期，即为ⅰb；
当t∈[85～140)，ds为中成岩阶段早期，即ⅱa；
当t∈[140～175)，ds为中成岩阶段晚期，即为ⅱb；
当t∈[175～200)，ds为晚成岩阶段，即为ⅲ；
3确定演化期次l3溶蚀作用
(1)将东营组二段l网格化，每个网格可用wi(x，y)表示；
(2)确定演化期次l3网格属性
a.根据步骤2确定网格wi(x，y)演化期次l3成岩阶段s_q；
b.根据研究区沉积相数据确定网格wi(x，y)演化期次l3的沉积相属性f_m；
c.根据研究区岩性数据确定网格wi(x，y)的岩性属性r_n；
d.根据研究区流体数据确定网格wi(x，y)演化期次l3的流体属性p_o；
e.根据上述步骤a、b、c和d确定每个网格wi(x，y)演化期次l3的综合属性g(s_q，f_m，r_n，p_o)，即：
wi(x，y)＝s_q+f_m+r_n+p_o；
(3)基于研究区地质资料，确定研究区已知井hj所在网格wj(x，y)不同演化期次的成岩阶段s_q_j以及成岩作用d_s_j；
(4)确定已知井hj所在网格wj(x，y)不同演化期次不同演化阶段的综合属性gj(f_m_j，r_n_j，p_o_j)，建立已知井hj所在网格wj(x，y)不同演化期次不同演化阶段成岩作用d_s与综合属性g(f_m_j，r_n_j，p_o_j)的对应关系，即为：
wj(s_q_j，d_s_j)＝f_m_j+r_n_j+p_o_j；
(5)任取一未知网格wi(x，y)，确定未知网格wi(x，y)成岩阶段s_q时的综合属性g(f_m，r_n，p_o)，即为：
gi＝(f_m_i，r_n_i，p_o_i)；
(6)将网格wi(x，y)成岩阶段s_q时的成岩作用综合属性gi与井hj所在网格wj(x，y)成岩阶段s_q的综合属性gj相比较，即为：
gi-gj＝(f_m_i，r_n_i，p_o_i)-(f_m_j+r_n_j+p_o_j)；
如果满足gi-gj＝0，则未知网格wi(x，y)的成岩作用d_s_i与井j具有相同的成岩作用d_s_j；
如果gi-gj≠0，则按不同属性优先级顺序，即为：沉积相f_m一级、岩性r_n二级、流体性质p_o三级对未知网格wi(x，y)的成岩作用进行判识，即：
ⅰ：f_m_i-f_m_j＝0，r_n_i-r_n_j≠0，p_o_i-p_o_j≠0；
ⅱ：f_m_i-f_m_j＝0，r_n_i-r_n_j＝0，p_o_i-p_o_j≠0；
未知网格wi(x，y)的成岩作用d_s_i满足条件i时，未知网格wi(x，y)与具有相同属性f_j的网格wj(x，y)的成岩作用相同；
d_s_i满足条件ⅱ时，未知网格wi(x，y)与具有相同属性f_m_j、r_n_j的网格wj(x，y)的成岩作用相同；
如在时间t＝5.0ma，埋深h＝250m处，网格w30＝(s_1，f_43，r_6，p_2)-(s_1_cfd17-3-1，f_m_cfd17-3-1，r_n_cfd17-3-1+p_o_cfd17-3-1)＝0，即网格w30与井cfd17-3-1具有相同的成岩作用，也即：w30＝d1-2_d4-2；
如在时间t＝3.6ma，埋深h＝180m处，网格w40＝(s_1，f_43，r_10，p_2)-(s_1_cfd18-2e-1，f_cfd18-2e-1，r_n_cfd18-2e-1，p_o_cfd18-2e-1)＝0，即网格w40与井cfd18-2e-1具有相同的溶蚀作用，也即：w40＝d1-2_d4-3；
如在时间t＝5.3ma，埋深h＝265m处，网格w50＝(s_1，f_26，r_11，p_2)网格属性满足条件ⅱ，即：f_m-f_m_bz-13-1-3＝0，r_n-r_n_bz-13-1-3＝0，p_o-p_o_bz-13-1-3≠0；则w50与井bz-13-1-3具有相同溶蚀作用，也即：w50＝d1-2_d6-3；
(7)根据每个网格wi(x，y)演化期次l3的成岩作用确定研究区成岩作用(图16)。
4)计算演化期次l3压实量、胶结量、溶蚀量
(1)根据成岩序列数据，确定已知井所在网格wj(x，y)的成岩阶段s_q所对应的压实作用、胶结作用、溶蚀作用下的压实级别数学模型a_k_j、胶结级别数学模型b_k_j、溶蚀级别数学模型c_k_j，即wj(s_q_j，d_s_j)＝wj(a_k_j，b_k_j，c_k_j)；
(2)将wi(x，y)的wi(s_q_j，d_s_j)与wj(s_q_j，d_s_j)比对，满足wi(s_q_j，d_s_j)＝wj(s_q_j，d_s_j)，则wi(x，y)与wj(s_q_j，d_s_j)具有相同的压实级别数学模型a_k、胶结级别数学模型b_k、溶蚀级别数学模型c_k，即a_k_i＝a_k_j、b_k_i＝b_k_j、c_k_i＝c_k_j；如果出现同一未知网格wi(x，y)与多个已知网格wj(s_q_j，d_s_j)相同，则a_k_i＝1/n∑a_k_j、b_k_i＝1/n∑b_k_j、c_k_i＝1/n∑c_k_j，n为相同网格wj(s_q_j，d_s_j)的网格数；
(3)根据网格wi(x，y)不同演化期次的压实级别数学模型a_k、胶结级别数学模型b_k、溶蚀级别数学模型c_k计算每个网格wi(x，y)不同演化期次的压实率rcom、胶结率rcem、溶蚀率rdis；
(4)根据每个网格wi(x，y)的压实率rcom、胶结率rmi_cem、溶蚀率rmi_dis计算压实量vcom，胶结量vmi_cem，溶蚀量vmi_dis，
如在时间t＝5.0ma，埋深h＝250m处，网格w30＝(s_1，f_43，r_6，p_2)-(s_1_cfd17-3-1，f_m_cfd17-3-1，r_n_cfd17-3-1+p_o_cfd17-3-1)＝0，即网格w30与井cfd17-3-1具有相同的压实级别数学模型a_k、胶结级别数学模型b_k，即a_k_i＝a_3_cfd17-3-1、b_k_i＝b_4_cfd17-3-1；
rcom＝4.3800*ln(dep)-1.9551＝4.3800*ln(250)-1.9551＝22.2
rcal_cem＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/5＝[-4.859e-06*(5)^4+0.001028*(5)^3-0.0745*(5)^2+2.862*(5)+10.6]/5＝4.6
φ＝φ0-vcom-vmi_cem+vmi_dis＝40-8.9-1.8＝29.3
如在时间t＝3.6ma，埋深h＝180m处，网格w40＝(s_1，f_43，r_10，p_2)-(s_1_cfd18-2e-1，f_cfd18-2e-1，r_n_cfd18-2e-1，p_o_cfd18-2e-1)＝0，即网格w40与井cfd18-2e-1具有相同的压实级别数学模型a_k、胶结级别数学模型b_k，即a_k_i＝a_3_cfd18-2e-1、b_k_i＝b_4_cfd18-2e-1；
rcom＝4.3800*ln(dep)-1.9551＝4.3800*ln(250)-1.9551＝20.8
rcal_cem＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/5＝[-4.859e-06*(3.6)^4+0.001028*(3.6)^3-0.0745*(3.6)^2+2.862*(3.6)+10.6]/5＝4.0
φ＝φ0-vcom-vmi_cem+vmi_dis＝40-8.3-1.6＝30.1
如在时间t＝5.3ma，埋深h＝265m处，网格w50＝(s_1，f_26，r_11，p_2)网格属性满足条件ⅱ，即：f_m-f_m_bz-13-1-3＝0，r_n-r_n_bz-13-1-3＝0，p_o-p_o_bz-13-1-3≠0；则w50与井bz-13-1-3具有相同压实级别数学模型a_3、胶结级别数学模型b_4、溶蚀级别数学模型c_k，即a_k_i＝a_2_bz-13-1-3、b_k_i＝b_2_bz-13-1-3；
rcom＝2.9200*ln(dep)-1.3034＝2.9200*ln(265)-1.3034＝15.0
rcal_cem＝[-4.859e-06*(x)^4+0.001028*(x)^3-0.0745*(x)^2+2.862*(x)+10.6]/8＝[-4.859e-06*(3.6)^4+0.001028*(3.6)^3-0.0745*(3.6)^2+2.862*(3.6)+10.6]/8＝2.5
φ＝φ0-vcom-vmi_cem+vmi_dis＝40-8.3-1.0＝30.7
(5)根据每个网格wi(x，y)演化期次l3的压实量、胶结量、溶蚀量，确定目标储层演化期次l3的孔隙分布(图19)；
5)确定演化期次l2、l1压实作用、胶结作用、溶蚀作用
利用步骤3，采用相同办法，确定东营组二段l演化期次l2和演化期次l1的压实作用、胶结作用、溶蚀作用(图17、图18)；
6)确定演化期次l2、l1压实量、胶结量、溶蚀量
利用步骤4，采用相同办法，确定东营组二段l演化期次l2和演化期次l1的压实量、胶结量、溶蚀量，确定目标储层演化期次l3和演化期次l1的孔隙度分布(图20、图21)；
表1研究区东营组二段压实作用表
表2研究区东营组二段胶结作用表
表3研究区东营组二段溶蚀作用表
其它未详细说明的部分均为现有技术。尽管上述实施例对本发明做出了详尽的描述，但它仅仅是本发明一部分实施例，而不是全部实施例，人们还可以根据本实施例在不经创造性前提下获得其他实施例，这些实施例都属于本发明保护范围。