本发明涉及油气勘探开发技术领域，具体地指一种碳酸盐岩储层成岩模拟孔隙度的定量预测方法。
背景技术：
碳酸盐储层在世界油气格局中占有相当大比例。从全球范围看，碳酸盐岩储层占沉积岩的20％，占油气探明储量的50％以上。在我国，近年来随着对海相碳酸盐岩储层油气勘探的不断深入，在渤海湾、塔里木、鄂尔多斯、四川盆地相继发现了多个大、中型油气田，展示了碳酸盐岩储层勘探的巨大潜力。由于碳酸盐岩储层沉积时代久远，经历了多次构造运动，遭受强烈的风化、剥蚀和淋滤作用，地层非均质性强，导致碳酸盐岩储层的定量评价成为碳酸岩盐油藏研究的主要技术难题之一。
目前，利用地震资料和测井资料，做波阻抗反演，拟合孔隙度与波阻抗的关系，构建地质数学模型，进而对储层孔隙度大小进行三维空间预测，取得了比较好的应用效果。但受地震资料分辨率的影响，此方法在薄层碳酸盐储层孔隙度预测上具有很大的局限性。
通过井上孔隙度大小分析测试数据，拟合孔隙度与深度之间的函数关系，再结合目的层埋深数据，对碳酸盐岩储层孔隙度大小空间分布作预测。但由于这种简单的数据拟合没有充分考虑不同层系之间沉积环境、岩性、成岩物质、成岩作用以及构造作用的不同，而只考虑了沉积环境、岩性、成岩物质以及成岩作用的综合结果，导致这种储层预测方法具有很大的局限性，不适合非均质性强的碳酸盐岩储层预测。
技术实现要素：
本发明针对现有技术的不足，提供了一种碳酸盐岩储层成岩模拟孔隙度的定量预测方法，该方法基于碳酸盐储层沉积相、成岩物质、沉积旋回、成岩作用以及地层厚度的研究，建立“地质参数-成岩作用-孔隙度”预测模型，进而对碳酸盐岩储层孔隙度大小进行横向预测，确定碳酸盐岩储层孔隙空间分布，为油藏评价提供依据。本研究充分考虑了碳酸盐岩储层非均质性的触发机制，能对碳酸盐岩储层的定量评价起到很好的预测作用。
为实现上述目的，本发明提供的一种碳酸盐岩储层成岩模拟孔隙度的定量预测方法，其特征在于：包括以下步骤：
1)收集研究区地质资料
2)原始资料数据库建立
建立的原始资料数据库包括以下几种：
(1)沉积相类型数据库f_m
(2)成岩相类型数据库d_s
(3)物源数据库o_n
(4)沉积旋回数据库c_p
(5)沉积厚度数据库t_h
3)建立孔隙度预测模型
(1)储层网格化
将研究区储层网格化；储层网格化是将空间上不均匀分布的数据按一定方法(如滑动平均法、克里格法或其他适当的数值推算方法)归算成规则网格中的代表值(趋势值)的过程。研究区储层的每个网格用wi(x，y，z)表示。
(2)确定网格属性；
a.根据研究区沉积相-岩相数据确定网格wi(x，y，z)的沉积相-岩相f_m；
b.根据研究区成岩相确定网格wi(x，y，z)的成岩相d_s；
c.根据研究区物源数据确定网格wi(x，y，z)的物源属性o_n；
d.根据研究区沉积旋回数据确定网格wi(x，y，z)的旋回属性c_p；
e.根据研究区构造等值线数据确定网格wi(x，y，z)的厚度属性t_h；
f.根据上述步骤a、b、c、d和e，确定每个网格的综合属性wi(x，y，z)＝g(f_m，d_s，o_n，c_p，t_h)；
(3)确定网格孔隙度预测模型；
a.确定相同沉积相-岩相区域f_m的孔隙度预测模型φ1i，
在研究区选择相同沉积相-岩相中的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到f_m区域的孔隙度预测模型，即为：
φ1i＝f1i(x，y，z，f_m)；
b.确定相同成岩相区域d_s的孔隙度预测模型φ2i，
在研究区选择相同成岩相区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到d_s区域的孔隙值预测模型，即为：
φ2i＝f2i(x，y，z，f_m)；
c.确定相同物源区域o_n的孔隙度预测模型φ3i，
在研究区选择相同物源区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到o_n区域的孔隙值预测模型，即为
φ3i＝f3i(x，y，z，o_n)；
d.确定相同沉积旋回区域c_p的孔隙度预测模型φ4i，
在研究区选择相同沉积旋回区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到c_p区域的孔隙值预测模型，即为：
φ4i＝f4i(x，y，z，c_p)；
e.确定相同厚度区域t_h的孔隙度预测模型φ5i，
在研究区选择相同厚度区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到t_h区域的孔隙值预测模型，即为
φ5i＝f5i(x，y，z，c_p)；
(4)孔隙模拟计算；
根据每个网格的综合属性wi(x，y，z)＝g(f_m，d_s，o_n，c_p，t_h)确定模拟点孔隙预测数学模型，并计算得到每个网格的孔隙度；
其中：φ为现今孔隙模拟值；φji为孔隙度预测模型，1≤j≤5；fj为孔隙度影响因子，fj＝0或fj＝1；kji为孔隙度影响系数。
进一步地，所述步骤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；
(2)成岩相类型数据库d_s
碳酸盐岩储层成岩相包括：胶结相d_1、白云化相d_2、方解石化相d_3、硅化相d_4、重结晶相d_5、强溶蚀相d_6、弱溶蚀相d_7、压溶相d_8、破裂相d_9、烃类充注相d_10；
(3)物源数据库o_n
物源数据库o_n包括物源东o_1、物源东西o_2、物源西o_3、物源西南o_4、物源南o_5、物源东南o_6、物源东o_7、物源东北o_8；
(4)沉积旋回数据库c_p
对碳酸盐岩储层沉积旋回类型进行概括和分类，其中包括：正旋回c_1、反旋回c_2、复合旋回c_3；
(5)沉积厚度数据库t_h
根据不同的沉积厚度将沉积厚度数据库t_h分为以下级别：
ⅰ级厚度t_1(10000≥t＞500米)、ⅱ级厚度t_2(500≤t＜200米)、ⅲ级厚度t_3(200≤t＜100米)、ⅳ级厚度t_4(100≤t＜50米)、ⅴ级厚度t_5(t≤50米)。
本发明的有益效果在于：
本发明基于碳酸盐储层沉积相、成岩物质、沉积旋回、地层厚度、成岩作用的研究，建立“地质参数-成岩作用-孔隙度”预测模型，对碳酸盐储层孔隙度进行横向预测，确定储层孔隙空间分布，由于该方法充分考虑了沉积环境、成岩物质、成岩作用以及构造作用的不同以及响应特征，对碳酸盐储层尤其是非均质性强的碳酸盐岩储层孔隙分布能起到很好预测作用。
附图说明
图1是研究区长兴组沉积相图；
图2是研究区长兴组成岩相图；
图3是研究区长兴组埋深等值线图；
图4是研究区长兴组物源分布图；
图5是研究区长兴组沉积旋回图；
图6是研究区长兴组地层厚度分区图；
图7是研究区长兴组孔隙分布图；
具体实施方式
为了更好地解释本发明，以下结合具体实施例进一步阐明本发明的主要内容，但本发明的内容不仅仅局限于以下实施例。
实施例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；
(2)成岩相类型数据库d_s
碳酸盐岩储层成岩相主要包括：胶结相d_1、白云化相d_2、方解石化相d_3、硅化相d_4、重结晶相d_5、强溶蚀相d_6、弱溶蚀相d_7、压溶相d_8、破裂相d_9、烃类充注相d_10。
(3)物源数据库o_n
物源数据库o_n包括物源东o_1、物源东西o_2、物源西o_3、物源西南o_4、物源南o_5、物源东南o_6、物源东o_7、物源东北o_8；
(4)沉积旋回数据库c_p
对碳酸盐岩储层沉积旋回类型进行概括和分类，其中包括：正旋回c_1、反旋回c_2、复合旋回c_3；
(5)沉积厚度数据库t_h
根据不同的沉积厚度将沉积厚度数据库t_h分为以下级别：
ⅰ级厚度t_1(t＞500米)、ⅱ级厚度t_2(500≤t＜200米)、ⅲ级厚度t_3(200≤t＜100米)、ⅳ级厚度t_4(100≤t＜50米)、ⅴ级厚度t_5(t≤50米)；
3)建立孔隙度预测模型
(1)将研究区储层网格化；储层网格化储层网格化是将空间上不均匀分布的数据按一定方法(如滑动平均法、克里格法或其他适当的数值推算方法)归算成规则网格中的代表值(趋势值)的过程。研究区储层的每个网格用wi(x，y，z)表示。
(2)确定网格属性；
a.根据研究区沉积相-岩相数据确定网格wi(x，y，z)的沉积相-岩相f_m；
b.根据研究区成岩相确定网格wi(x，y，z)的成岩相d_s；
c.根据研究区物源数据确定网格wi(x，y，z)的物源属性o_n；
d.根据研究区沉积旋回数据确定网格wi(x，y，z)的旋回属性c_p；
e.根据研究区构造等值线数据确定网格wi(x，y，z)的厚度属性t_h；
f.根据上述步骤a、b、c、d和e，确定每个网格的综合属性wi(x，y，z)＝g(f_m，d_s，o_n，c_p，t_h)；
(3)确定网格孔隙度预测模型；
a.确定相同沉积相-岩相区域f_m的孔隙度预测模型φ1i，
在研究区选择相同沉积相-岩相中的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到f_m区域的孔隙度预测模型，即为：
φ1i＝f1i(x，y，z，f_m)；
b.确定相同成岩相区域d_s的孔隙度预测模型φ2i，
在研究区选择相同成岩相区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到d_s区域的孔隙值预测模型，即为：
φ2i＝f2i(x，y，z，f_m)；
c.确定相同物源区域o_n的孔隙度预测模型φ3i，
在研究区选择相同物源区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到o_n区域的孔隙值预测模型，即为
φ3i＝f3i(x，y，z，o_n)；
d.确定相同沉积旋回区域c_p的孔隙度预测模型φ4i，
在研究区选择相同沉积旋回区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到c_p区域的孔隙值预测模型，即为：
φ4i＝f4i(x，y，z，c_p)；
e.确定相同厚度区域t_h的孔隙度预测模型φ5i，
在研究区选择相同厚度区域的单井分析测试数据，拟合孔隙值φ与深度z的函数关系，得到t_h区域的孔隙值预测模型，即为
φ5i＝f5i(x，y，z，c_p)；
(4)孔隙模拟计算；
根据每个网格的综合属性wi(x，y，z)＝g(f_m，d_s，o_n，c_p，t_h)确定模拟点孔隙预测数学模型，并计算得到每个网格的孔隙度；
其中：φ为现今孔隙模拟值；φji为孔隙度预测模型，1≤j≤5；fj为孔隙度影响因子，fj＝0或fj＝1；kji为孔隙度影响系数。
实施例2
基于上述方法以某省盆地东北部长兴组碳酸盐岩储层为研究对象，预测长兴组储层孔隙度空间分布，从而分析长兴组储层孔隙发育情况，具体如下：
1)收集研究区地质资料，包括以下内容：
(1)研究区长兴组沉积相图(图1)；
(2)研究区长兴组成岩相图(图2)；
(3)研究区长兴组埋深等值线图(图3)；
(4)研究区长兴组物源分布图(图4)；
(5)研究区长兴组沉积旋回图(图5)；
(6)研究区长兴组地层厚度分区图(图6)；
(7)研究区长兴组孔隙分布图(图7)；
(8)研究区单井孔隙度数据(表1、表2、表3)；
表1
表2
表3
2)建立孔隙度预测模型
(1)储层网格化
将研究区长兴组碳酸盐岩储层进行网格化，每个网格可用wi(x，y，z)表示。
(2)确定网格属性；
a.根据研究区沉积相数据，确定网格wi(x，y，z)的沉积相f_m；
b.根据研究区成岩相，确定网格wi(x，y，z)的成岩相d_s；
c.根据研究区物源数据，确定网格wi(x，y，z)的物源属性o_n；
d.根据研究区沉积旋回数据，确定网格wi(x，y，z)的旋回属性c_p；
e.根据研究区构造等值线数据，确定网格wi(x，y，z)的厚度属性t_h；
f.根据a、b、c、d、e确定每个网格的综合属性wi(x，y，z)＝g(f_m，d_s，o_n，c_p，t_h)；
(3)确定网格孔隙度预测模型；
a.确定相同沉积相区域f_m的孔隙度预测模型φ1i，即φ1i＝f1i(x，y，z，f_m)
根据井c16-1、c18-1、c25-1、c18-3、c23-2、c4-1、c17-1、c1-4、c18-2、c11-3、c21-1、c1-3、c8-3的孔隙度分析测试数据，确定开阔台地的孔隙度预测模型：
φ11＝0.0000000002688*h^4-0.000005183*h^3+0.03737*h^2-119.4*h+142640
其中，φ11为开阔台地相孔隙度，h为埋深；
根据井c16-3、c5-1、c23-1、c2-1、c8-2确定台地边缘的孔隙度预测模型：
φ12＝-0.00000000000303*h^4+0.00000005376*h^3-0.0003574*h^2+1.056*h-1154.5
其中，φ12为台地边缘相孔隙度，h为埋深；
根据井c16-2、c4-1、c9-1、c9-2、c1-2、c1-1、c3-1确定台地斜坡的孔隙度预测模型：
φ13＝-0.000000000000866*h^4+0.00000001583*h^3-0.0001083*h^2+0.3288*h-360
其中，φ13为台地斜坡相孔隙度，h为埋深；
b.确定相同成岩相区域d_s的孔隙度预测模型φ2i，即φ2i＝f2i(x，y，z，d_s)
根据井c16-1、c17-1、c18-2、c11-3、c3-1、c1-2、c9-2、c25-1、c18-3、c4-1、c2-1、c23-1、c1-4、c8-2、c7-2的孔隙度分析测试数据，确定溶蚀作用的孔隙度预测模型：
φ21＝0.000000000006707*h^4-0.0000001308*h^3+0.0009544*h^2-3.091*h+3765
其中，φ21为溶蚀相孔隙度，h为埋深；
根据井c16-3、c5-1确定胶结作用的孔隙度预测模型：
φ22＝-0.0000000000007944*h^4+0.00000001329*h^3-0.00008324*h^2+0.2315*h-228.4576
其中，φ22为胶结相孔隙度，h为埋深；
根据井c1-3、c8-3、c1-1确定压溶作用相的孔隙度预测模型：
φ23＝0.00000000001856*h^4-0.0000003715*h^3+0.002782*h^2-9.243*h+11516
其中，φ23为压溶相孔隙度，h为埋深；
根据井c18-1、c23-2、c16-2、c21-1的孔隙度分析测试数据，确定方解石化的孔隙度预测模型：
φ24＝-0.00000000000001261*h^4+0.0000000002355*h^3-0.0000001649*h
^2+0.005129*h-25.934
其中，φ24为方解石化相孔隙度，h为埋深；
c.确定相同物源区域o_n的孔隙度预测模型φ3i，即φ3i＝f3i(x，y，z，o_n)
根据井c16-1、c18-1、c17-1、c23-2、c25-1、c18-3、c16-3的孔隙度分析测试数据，确定北西物源的孔隙度预测模型：
φ31＝-0.0000000000007794*h^4+0.0000000131*h^3-0.00008183*h^2+0.2249*h-211
其中，φ31为北西物源孔隙度，h为埋深；
根据井c2-1、c5-1、c4-1、c23-1、c16-2、c1-4、c9-1、c18-2、c11-3、c21-1、c9-2、c1-3、c8-3、c3-1、c8-2、c1-1的孔隙度分析测试数据，确定东南物源的孔隙度预测模型：
φ32＝-0.000000000005794*h^4+0.0000001048*h^3-0.0007093*h^2+2.13*h-2372
其中，φ32为东南物源孔隙度，h为埋深；
d.确定相同沉积旋回区域c_p的孔隙度预测模型φ4i，即φ4i＝f4i(x，y，z，c_p)
根据井c16-1、c18-1、c25-1、c16-3、c18-3、c16-2、c23-2、c17-1的孔隙度分析测试数据，确定正旋回的孔隙度预测模型：
φ41＝0.00000000000358*h^4-0.00000006975*h^3+0.0005085*h^2-1.644*h+2000
其中，φ41为正旋回孔隙度，h为埋深；
根据井c2-1、c5-1、c4-1、c23-1、c1-4、c9-1、c18-2、c11-3的孔隙度分析测试数据，确定复合旋回的孔隙度预测模型：
φ42＝-0.0000000000003143*h^4+0.000000005504*h^3-0.00003613*h^2+0.1054*h-100.4
其中，φ42为复合旋回孔隙度，h为埋深；
根据井c3-1、c1-2、c21-1、c9-2、c7-2、c1-1、c8-3、c1-3、c8-2的孔隙度分析测试数据，确定反旋回的孔隙度预测模型：
φ43＝0.00000000001415*h^4-0.0000002757*h^3+0.002013*h^2-6.529*h+7948.8
其中，φ43为反旋回孔隙度，h为埋深；
e.确定相同厚度区域t_h的孔隙度预测模型φ5i，即φ5i＝f5i(x，y，z，c_p)
根据井c16-1、c17-1、c23-2的孔隙度分析测试数据，确定ⅱ级厚度孔隙度预测模型：
φ51＝0.000000000003449*h^4-0.00000006874*h^3+0.0005132*h^2-1.7*h+2121
其中，φ51为ⅱ级厚度孔隙度，h为埋深；
根据井c18-1、c25-1、c18-3、c16-2、c18-2、c1-4、c9-1、c9-2的孔隙度分析测试数据，确定ⅲ级厚度孔隙度预测模型：
φ52＝-0.0000000000006364*h^4+0.00000001141*h^3-0.00007664*h^2+0.2285*h-241.6
其中，φ52为ⅲ级厚度孔隙度，h为埋深；
根据井c16-3、c2-1、c5-1、c4-1、c23-1、c21-1、c11-3、c1-2、c1-1、c3-1、c7-2、c8-3、c1-3、c8-2的孔隙度分析测试数据，确定ⅳ级厚度孔隙度预测模型：
φ53＝-0.00000000000141*h^4+0.0000000251*h^3-0.0001673*h^2+0.4949*h-530
其中，φ53为ⅳ级厚度孔隙度，h为埋深；
(4)孔隙模拟计算；
a.根据研究区长兴组碳酸盐岩储层储层属性综合分布图，确定研究区储层每个网格的综合属性wi(x，y，z)＝g(f_m，d_s，o_n，c_p，t_h)；
b.根据步骤(2)中a、b、c、d、e确定每个模拟点孔隙预测数学模型，计算其孔隙度，即：
其中：φ为现今孔隙模拟值；φji为孔隙度预测模型，1≤j≤5；fj为孔隙度影响因子，fj＝0或fj＝1，在本实例中fj＝1；kji为孔隙度影响系数，在本实例中k1i＝0.4、k2i＝0.3、k3i＝0.2、k4i＝0.15、k5i＝0.3；
如网格w10，沉积相为开阔台地相，成岩相为溶蚀相，物源方向北西，沉积旋回为正旋回，沉积厚度属于ⅱ级厚度，埋深4800米，即网格w10＝(f_16，d_6，o_2，c_1，t_2)；网格w10的孔隙度预测模型如下：
如网格w20，沉积相为台地边缘相，成岩相为溶蚀相，物源方向东南，沉积旋回为复合，沉积厚度属于ⅲ级厚度，埋深5000米，即网格w20＝(f_11，d_6，o_6，c_3，t_3)；网格w20的孔隙度预测模型如下：
如网格w30，沉积相为台地边缘相，成岩相为压溶相，物源方向东南，沉积旋回为复合，沉积厚度属于ⅳ级厚度，埋深5100米，即网格w30＝(f_16，d_8，o_6，c_3，t_4)；网格w30的孔隙度预测模型如下：
如网格w40，沉积相为台地斜坡相，成岩相为方解石化，物源方向北西，沉积旋回反旋回，沉积厚度属于ⅲ级厚度，埋深5250米，即网格w40＝(f_9，d_3，o_2，c_2，t_3)；网格w40的孔隙度预测模型如下：
c.根据每个空间网格的孔隙度，得到研究区长兴组碳酸盐岩储层孔隙度分布(图7)。
其它未详细说明的部分均为现有技术。尽管上述实施例对本发明做出了详尽的描述，但它仅仅是本发明一部分实施例，而不是全部实施例，人们还可以根据本实施例在不经创造性前提下获得其他实施例，这些实施例都属于本发明保护范围。