[发明专利]一种基于数字岩心模型的储层敏感性伤害模拟方法

权利要求书

1.一种基于数字岩心模型的储层敏感性伤害模拟方法，其特征在于，包括以下步骤：

步骤1，基于真实储层二维信息，利用改进混合算法和聚类算法构建含多种岩石矿物的数字岩心模型；

步骤2，通过室内实验得到不同模拟条件下，不同类型岩石矿物的体积变化、转化规律的实验参数；

步骤3，基于不同模拟条件下各类岩石矿物的体积膨胀系数，计算得到不同类型岩石矿物基团的膨胀体积，以步骤1中的含多种岩石矿物组分的数字岩心模型为基础，结合不同类型岩石矿物的属性及形态学中的膨胀算法，对不同模拟条件下岩石矿物的膨胀过程进行模拟；

步骤4，基于不同模拟条件下各类岩石矿物的体积溶蚀系数，计算得到不同类型岩石矿物基团的溶蚀体积，以步骤1中的含多种岩石矿物组分的数字岩心模型为基础，结合不同类型岩石矿物的属性及形态学中的腐蚀算法，对不同模拟条件下岩石矿物的溶蚀过程进行模拟；

步骤5，基于不同模拟条件下各类岩石矿物的转化率，计算得到不同类型岩石矿物基团相互转化的体积，以步骤1中的含多种岩石矿物组分的数字岩心模型为基础，优先考虑岩石矿物基团较小的岩石矿物基团进行转化，包括转化相岩石矿物的减少和被转化相岩石矿物增加两个过程，对不同条件下岩石矿物的转化过程进行模拟；

其中：

所述步骤1中，真实储层二维信息包括铸体薄片、岩石粒度分布、黏土矿物分布、黏土矿物产状特征；具体的模型构建步骤包括：

第一步，利用过程法构建基础数字岩心模型时，考虑黏土矿物的总含量，在沉积过程中，根据真实储层的粒度分布情况，随机选择沉积颗粒的半径，沉积颗粒的尺寸不但由原始的沉积颗粒粒度分布决定，同时额外考虑黏土矿物与储层砂岩颗粒之间的比例，在满足高能环境和重力势能梯度最大的下落模拟原则的基础上模拟沉积过程，并结合真实岩心孔隙度，选择压实因子控制数字岩心模型的孔隙度；

第二步，将单位体像素点的空间占位，即点、线和面占位三种类型，按其对邻域不稳定性的贡献程度赋予权值，其中面为5，边为3，点为2；在选取交换单位体像素点时，计算该体像素点与邻域占位点、线和面上的不稳定性贡献程度S，并基于模拟退火算法中能量值下降的过程，引入交换单位体像素点对其邻域不稳定性的贡献程度参数S_d，对交换点的可交换性进行判断，提高交换单位体像素点的有效性，其中S_d为与模拟过程中系统能量相关的无因次值：

S_d＝N×β(E₀-E_i/ΔE_max) (1)

式中，N为单位体像素点影响的邻域接触点的个数，无量纲；β为单位体像素点对邻域不稳定性系数，无量纲；E₀为系统的初始能量，无量纲；E_i为第i次降温后系统的能量，无量纲；ΔE_max为初始模型和基于储层岩石二维信息的参考模型系统的能量差值，无量纲，初始模型是指过程法构建的基础数字岩心模型；

第三步，利用改进混合算法构建初始数字岩心模型的步骤为：

①建立基于储层岩石二维信息的参考模型，将过程法构建的基础数字岩心模型作为改进混合算法的初始模型，设定初始温度，并计算初始系统的相关参数，包含自相关函数、线性路径函数、分形特征函数和能量值；

②在保证模拟退火降温过程随机性的基础上，计算交换单位体像素点26个空间占位对邻域不稳定性的贡献程度S；当SS_d时，认为该点的不稳定程度较高，可作为系统更新的交换点；当SS_d时，则重复步骤②；

③计算交换单位体像素点后系统的相关参数，包括单点概率函数、自相关函数、线性路径函数、分形函数和能量值，计算与未交换前系统的能量差值ΔE；当ΔE0时，更新系统；当ΔE0时，根据Metropolis准则来判断系统是否更新，即在一定的概率条件下接受系统更新；如果判断后不满足系统更新条件，则返回步骤②；

④判断内循环终止条件，即判断在同一温度条件下系统能量差值是否小于设定最小能量差值；同时为避免系统刚降温，系统能量上升而立刻导致内循环结束而产生的降温，通过设定系统更新的失败率f_f来避免该现象的出现，其中：

式中，N_f为导致系统能量回升的更新失败的次数；N为系统更新的总次数；

当f_f大于一定值后，则进行降温处理，降温过程采取等比降温方案，并返回步骤②；

⑤当模拟过程温度降低到最终设定温度时或与上次降温的系统能量差值ΔE小于设定值时，整个模拟过程终止；

作为约束条件，模拟退火算法中使用的统计函数包括：单点概率函数P(r)、自相关函数、线性路径函数和分形函数，利用自相关函数和线性路径函数对初始系统进行退火模拟，当模型具备一定分形特征后，引入分形函数进一步约束重建模型；

第四步，将混合算法重建后初始数字岩心模型中的类球岩石颗粒，与过程法中构建的基础数字岩心模型的原始球形岩石颗粒相比较并取二者补集，将初始数字岩心模型初步划分为岩石骨架相、孔隙相和黏土矿物相三大类；

第五步，通过Hoshen-Kopelman算法对初始数字岩心模型中的黏土矿物基团进行统计和划分，其中被M相占据的概率为c，被T相占据的概率为1-c，对于晶格中的每一个占位i，当其被M相占据时，则给该占位赋予一个基团标记其中α是基团标记的特征符号，t为基团标记的次数，某一离散点的标记由一系列自然数表示：

在这一系列自然数中只有一个自然数是基团α的准确标记，该标记为且该值是集合(3)中所有自然数的最小值，其它各基团标记之间的关系则由以下整数集给出：

其中，只有是正整数元素，该值为基团中M相的个数，当进行第t次标记时，若基团中M相个数少于上次标记过程基团α的M相个数，则将该差值表示为相应t次的基团α的T相个数，(4)中的其它元素皆为负整数，反映了与其它基团标记的关系，与的关系用式(5)表示：

检查被判断离散点是否有被扫描过的相邻离散点，若相邻离散点为T相，则将当前被判断离散点赋予新基团的标记；如果有一个相邻离散点已经赋予基团标记，则将当前网格与相邻离散点赋予相同的标记；如果有一个以上的相邻离散点已经赋予基团标记，且基团标记各不相同，则将基团中所有离散点赋予相同的标记，最后统计并划分模型中黏土矿物相基团的个数及尺寸；

第六步，较大的连通基团为黏土相中基团尺寸大于相邻基质颗粒尺寸的黏土矿物基团，通过K-means算法对初始数字岩心模型中黏土矿物相基团尺寸较大的黏土矿物基团进行划分，具体步骤如下：

①读取数据样本的集合；

②设定样本聚类的个数k，随机的选取k个数据样本作为初始的数据样本聚类中心；

③计算欧氏距离，计算数据样本中每个数据到各聚类中心的欧式几何距离，然后根据最小误差平方和准则函数将数据按照远近距离划分到相应的不同聚类中心所对应的聚类当中；

④更新聚类中心，将每个聚类中所有数据的均值作为各个聚类新的中心，并以最小误差平方和准则重新计算新的聚类中心的值；

⑤迭代判别，将步骤④中计算得到的数值与前一次计算得到的数值相比较，如果两者差值小于或等于预先设定的临界值，则停止迭代，否则重新进行步骤③进行迭代；

⑥输出数据样本及聚类结果，包括每个聚类的聚类中心、大小；

第七步，当黏土矿物基团边界的离散点为单个岩石颗粒时，则将该黏土矿物基团划分为交代形式，交代形式主要分布于岩石颗粒内，呈单个离散点的形式分布；当黏土矿物基团边界的相邻离散点为单个岩石骨架颗粒及孔隙时，则将该黏土矿物相基团划分为颗粒表面充填形式；

当黏土矿物基团边界的相邻离散点为多个岩石骨架颗粒及孔隙时，则将该黏土矿物基团划分为粒间充填形式；

将交代形式、颗粒表面充填形式和粒间充填形式的黏土矿物基团分别标记为A、B、C；最终得到不同结构黏土矿物基团分布和不同类型的黏土矿物基团分布；

第八步，基于Hoshen-Kopelman算法和K-means算法得到初始数字岩心模型中黏土矿物基团大小及数量分布，以及按结构划分得到的黏土矿物基团类型及数量分布，结合真实储层黏土含量及分布以及主要的黏土矿物结构特点，按黏土矿物相基团大小和结构特点将模型中的黏土矿物赋予相应的黏土性质，得到含多组分岩石矿物分布的数字岩心模型；

所述步骤2中，室内实验包括为模拟原始储层条件及不同生产条件下，当储层温度、压力、层内流体和注入流体性质发生变化时，不同岩石矿物的转化及体积变化情况，具体包括不同类型岩石矿物在不同模拟条件即不同温度，不同润湿环境条件，不同水类型下的膨胀率、溶蚀率和转化率。