本发明涉及地质及石油勘探技术领域，尤其涉及一种基于多点地质统计学的地质随机反演方法及装置。
背景技术：
目前，储层物性参数(包括孔隙度、泥质含量、含水饱和度等)是进行储层预测与评价、估算油气储量和确定开发井位的重要依据。地震数据中包含十分丰富的储层物性信息，具有横向分辨率较高的特征，是储层表征和预测中使用的重要资料，但是由于储层物性参数与地震数据之间的关系异常复杂，为利用地震数据获得储层物性参数带来巨大的困难。地震反演技术在储层表征和建模中能够发挥重要的作用。通常，地震反演可以分为确定性反演和随机反演。确定性反演由于地震资料的带限性质，影响了分辨率的进一步提高；而随机反演则与确定性反演相比，以地质统计为基础，测井数据为基本条件数据，在分辨率提高和不确定性降低上具有一定的优势。
地质统计学是随机反演的一种重要原理。两点地质统计学通过变差函数表征区域化变量的空间相关性，较为广泛的应用于科学研究的众多领域。在地球物理领域，序贯高斯模拟和序贯指示模拟是两点地质统计学中最常用的两种方法，主要用来处理连续变量(如：孔隙度、泥质含量等)和离散变量(如：岩相)的建模等问题。两点地质统计学方法计算效率快，能对离散变量和连续变量进行模拟和反演，但两点地质统计学方法利用变差函数描述地质变量的相关性和变异性，仅仅通过某个方向上两点之间的地质变量的变化关系来描述空间的变化特性，对目标体形态的再现能力较差。而多点地质统计学模拟作为两点地质统计学的补充和发展应运而生。多点地质统计学方法以像元为基本模拟单元，同时结合了基于目标的算法的优点，利用“训练图像”代替变差函数描述区域化变量的空间相关性，能够更好的刻画目标体的几何形态。目前，由于连续变量的训练图像难以获取，多点地质统计学主要用来获得离散变量(如岩相)的空间模型。多点地质统计学随机反演通常的做法是先在多点地质统计学模拟阶段，产生大量的等可能的储层岩相模型。根据这些岩相模型，从对应的储层物性参数的累积概率分布中抽样获得多个伪物性参数分布；接下来利用岩石物理理论将其转换为弹性信息；然后，正演合成地震记录与实际地震记录做匹配，若不能满足要求，则重新进行模拟一般需要模拟成百次才能获得比较满意的反演结果，计算量非常大，同时需要进行大量的正演模拟，巨大的计算量导致计算效率十分低效，这对多点地质统计学在储层表征和建模中的应用造成了阻碍。
可见，目前基于多点地质统计学的随机反演方法存在适用范围局限、计算效率低以及仅能获得岩相反演结果等问题。
技术实现要素：
本发明的实施例提供一种基于多点地质统计学的地质随机反演方法及装置，以解决目前基于多点地质统计学的随机反演方法存在适用范围局限、计算效率低以及仅能获得岩相反演结果的问题。
为达到上述目的，本发明采用如下技术方案：
一种基于多点地质统计学的地质随机反演方法，包括：
步骤1、通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数；
步骤2、获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型；
步骤3、统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型；
步骤4、根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型；
步骤5、对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录；
步骤6、将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果；
步骤7、将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
具体的，所述通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数，包括：
步骤101、建立研究区域的地质的模拟网格，获得训练图像，将研究区域的地质的条件数据分配至训练图像的网格节点上，并设置一随机路径；
步骤102、根据预先设置的数据模板扫描所述训练图像获得数据事件；所述数据模板为网格节点之间的空间组合关系；
步骤103、根据所述数据事件在所述训练图像中检索，得到待估计点处的条件概率分布；所述待估计点的条件概率分布为：
P{F(u)＝Fk|dn}≈ck(dn)/c(dn)
其中，c(dn)为所述数据事件在训练图像中的有效训练图像中重复的数量；ck(dn)为数据事件的所有重复的数量中F(u)＝Fk的个数；P{F(u)＝Fk|dn}为待估计点u的条件概率分布；
步骤104、根据所述待估计点的条件概率分布进行抽样，获得一模拟节点的模拟岩相，作为一次模拟结果，并将所述模拟结果加入到所述条件数据中；
步骤105、随机选取待估计点中的另一模拟节点，重复步骤103至步骤105，直至各待估计点均模拟结束，获得所述研究区域的初始岩相模型；
步骤106、获取预先设置的最大模拟次数。
具体的，所述获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型，包括：
步骤201、获取预先设置的最大扰动次数和权重因子；
步骤202、对所述初始岩相模型进行概率扰动更新，确定概率扰动更新后的联合概率；
所述概率扰动更新后的联合概率为：P(Fk|dobs)＝(1-r)*i0(u,fk)+r*P(Fk)；其中，为指示变量；P(Fk)表示空间中某未知点u处的岩相分布先验概率；P(Fk|dobs)为该未知点u的概率扰动更新后的联合概率；r表示权重因子，所述权重因子为采用非线性变换的权重因子；
步骤203、将所述概率扰动更新后的联合概率P(Fk|dobs)和预先得到的多点概率P(Fk|G)融合，确定扰动融合后的条件概率分布；
所述扰动融合后的条件概率分布为
其中τ1和τ2为模型参数，a＝(1-P(Fk))/P(Fk)，b＝(1-P(Fk|G))/P(Fk|G)，c＝(1-P(Fk|dobs))/P(Fk|dobs)；
步骤204、将扰动融合后的条件概率分布转化为累积概率分布，根据所述累积概率分布，采用随机抽样方法依次对模拟节点进行抽样，生成概率扰动更新后的岩相模型。
具体的，所述统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型，包括：
步骤301、根据预先获知的测井数据，统计不同岩相中的储层物性参数的分布特征，获得不同岩相中的储层物性参数的分布函数；
步骤302、根据随机抽样方法对概率扰动更新后的岩相模型中的每个网格节点进行随机抽样，获得每个相应网格节点处的物性参数状态值，由每个相应网格节点处的物性参数状态值构成所述伪储层物性参数模型。
具体的，所述根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型，包括：
根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系：e＝fRPM(r)+ε，将所述伪储层物性参数模型转换为弹性参数模型；
其中，fRPM(r)为通过理论岩石物理模型或者经验岩石物理关系获得的确定性岩石物理模型；e为弹性参数；r为储层物性参数；ε为确定性岩石物理模型的误差项。
具体的，所述对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录，包括：
步骤501、对所述弹性参数模型进行时深转换获得时间域弹性参数模型；
步骤502、根据所述时间域弹性参数模型计算反射系数；
步骤503、采用褶积模型正演合成地震记录其中，w为采用褶积模型正演合成的地震记录；b为从实际地震记录中提取的用于合成地震记录的地震子波；R为所述反射系数。
具体的，所述将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果，包括：
步骤601、将所述地震记录与预先确定的实际地震记录进行对比，确定所述地震记录与预先确定的实际地震记录的误差E＝||dobs-dsyn||2；其中，dobs为所述实际地震记录；dsyn为正演合成的地震记录；
步骤602、判断所述误差E是否小于或等于预先设置的误差容忍度；
步骤603、在所述误差E小于或等于所述误差容忍度时，将概率扰动更新后的岩相模型作为所述岩相模型结果，将所述伪储层物性参数模型作为伪储层物性参数模型结果；
步骤604、在所述误差E大于所述误差容忍度时，判断进行概率扰动更新的次数是否达到所述最大扰动次数；
步骤605、在进行概率扰动更新的次数未达到所述最大扰动次数时，返回步骤2的进行概率扰动更新的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度；
步骤606、在进行概率扰动更新的次数达到所述最大扰动次数时，判断当前进行模拟的次数是否达到所述最大模拟次数；
步骤607、在当前进行模拟的次数未达到所述最大模拟次数时，修改随机路径，并返回步骤1的通过多点地质统计学随机模拟算法确定初始岩相模型的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度；
步骤608、在当前进行模拟的次数达到所述最大模拟次数时，获取与预先设置的误差容忍度最接近的误差E对应的概率扰动更新后的岩相模型作为所述岩相模型结果，获取与预先设置的误差容忍度最接近的误差E对应的伪储层物性参数模型作为伪储层物性参数模型结果。
具体的，所述将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果，包括：
步骤701、获取待处理的伪储层物性参数模型结果，设置一初始体系动能CΓ(t)和退火终止条件；
步骤702、根据公式计算目标函数E(m(k))；其中，dj为观测到的地震数据；m(k)为第k次迭代的模型参数；
步骤703、根据公式m(k+1)＝m(k)+ξα对所述模型参数进行更新，确定第k+1次迭代的模型参数m(k+1)；ξ为0-1之间的随机数，α为步长；
步骤704、根据公式ΔE＝E(m(k+1))-E(m(k))计算体系势能ΔE；
步骤705、在ΔE<0时，采用m(k+1)替换m(k)；
步骤706、在ΔE≥0时，根据体系温度为T时的概率接受准则ρ(H)＝(e-H/PT)P，在m(k+1)和m(k)中确定概率ρ(H)的较大值对应的模型参数；其中，H＝ΔE+CΓ(t)；P为模型参数总数；
步骤707、在当前迭代的模型参数满足退火终止条件时，将m(k+1)作为储层物性参数反演结果；
步骤708、在当前迭代的模型参数不满足所述退火终止条件时，返回步骤703的模型参数更新过程，直至迭代的模型参数满足所述退火终止条件。
一种基于多点地质统计学的地质随机反演装置，包括：
初始岩相模型确定单元，用于通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数；
概率扰动更新单元，用于获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型；
伪储层物性参数模型转化单元，用于统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型；
弹性参数模型转换单元，用于根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型；
正演合成单元，用于对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录；
结果确定单元，用于将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果；
退火优化单元，用于将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
此外，所述初始岩相模型确定单元，具体用于执行：
步骤101、建立研究区域的地质的模拟网格，获得训练图像，将研究区域的地质的条件数据分配至训练图像的网格节点上，并设置一随机路径；
步骤102、根据预先设置的数据模板扫描所述训练图像获得数据事件；所述数据模板为网格节点之间的空间组合关系；
步骤103、根据所述数据事件在所述训练图像中检索，得到待估计点处的条件概率分布；所述待估计点的条件概率分布为：
P{F(u)＝Fk|dn}≈ck(dn)/c(dn)
其中，c(dn)为所述数据事件在训练图像中的有效训练图像中重复的数量；ck(dn)为数据事件的所有重复的数量中F(u)＝Fk的个数；P{F(u)＝Fk|dn}为待估计点u的条件概率分布；
步骤104、根据所述待估计点的条件概率分布进行抽样，获得一模拟节点的模拟岩相，作为一次模拟结果，并将所述模拟结果加入到所述条件数据中；
步骤105、随机选取待估计点中的另一模拟节点，重复步骤103至步骤105，直至各待估计点均模拟结束，获得所述研究区域的初始岩相模型；
步骤106、获取预先设置的最大模拟次数。
此外，所述概率扰动更新单元，具体用于执行：
步骤201、获取预先设置的最大扰动次数和权重因子；
步骤202、对所述初始岩相模型进行概率扰动更新，确定概率扰动更新后的联合概率；
所述概率扰动更新后的联合概率为：P(Fk|dobs)＝(1-r)*i0(u,fk)+r*P(Fk)；其中，为指示变量；P(Fk)表示空间中某未知点u处的岩相分布先验概率；P(Fk|dobs)为该未知点u的概率扰动更新后的联合概率；r表示权重因子，所述权重因子为采用非线性变换的权重因子；
步骤203、将所述概率扰动更新后的联合概率P(Fk|dobs)和预先得到的多点概率P(Fk|G)融合，确定扰动融合后的条件概率分布；
所述扰动融合后的条件概率分布为
其中τ1和τ2为模型参数，a＝(1-P(Fk))/P(Fk)，b＝(1-P(Fk|G))/P(Fk|G)，c＝(1-P(Fk|dobs))/P(Fk|dobs)；
步骤204、将扰动融合后的条件概率分布转化为累积概率分布，根据所述累积概率分布，采用随机抽样方法依次对模拟节点进行抽样，生成概率扰动更新后的岩相模型。
此外，所述伪储层物性参数模型转化单元，具体用于执行：
步骤301、根据预先获知的测井数据，统计不同岩相中的储层物性参数的分布特征，获得不同岩相中的储层物性参数的分布函数；
步骤302、根据随机抽样方法对概率扰动更新后的岩相模型中的每个网格节点进行随机抽样，获得每个相应网格节点处的物性参数状态值，由每个相应网格节点处的物性参数状态值构成所述伪储层物性参数模型。
此外，所述弹性参数模型转换单元，具体用于：
根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系：e＝fRPM(r)+ε，将所述伪储层物性参数模型转换为弹性参数模型；
其中，fRPM(r)为通过理论岩石物理模型或者经验岩石物理关系获得的确定性岩石物理模型；e为弹性参数；r为储层物性参数；ε为确定性岩石物理模型的误差项。
此外，所述正演合成单元，具体用于执行：
步骤501、对所述弹性参数模型进行时深转换获得时间域弹性参数模型；
步骤502、根据所述时间域弹性参数模型计算反射系数；
步骤503、采用褶积模型正演合成地震记录其中，w为采用褶积模型正演合成的地震记录；b为从实际地震记录中提取的用于合成地震记录的地震子波；R为所述反射系数。
此外，所述结果确定单元，具体用于执行：
步骤601、将所述地震记录与预先确定的实际地震记录进行对比，确定所述地震记录与预先确定的实际地震记录的误差E＝||dobs-dsyn||2；其中，dobs为所述实际地震记录；dsyn为正演合成的地震记录；
步骤602、判断所述误差E是否小于或等于预先设置的误差容忍度；
步骤603、在所述误差E小于或等于所述误差容忍度时，将概率扰动更新后的岩相模型作为所述岩相模型结果，将所述伪储层物性参数模型作为伪储层物性参数模型结果；
步骤604、在所述误差E大于所述误差容忍度时，判断进行概率扰动更新的次数是否达到所述最大扰动次数；
步骤605、在进行概率扰动更新的次数未达到所述最大扰动次数时，返回概率扰动更新单元的进行概率扰动更新的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度；
步骤606、在进行概率扰动更新的次数达到所述最大扰动次数时，判断当前进行模拟的次数是否达到所述最大模拟次数；
步骤607、在当前进行模拟的次数未达到所述最大模拟次数时，修改随机路径，并返回初始岩相模型确定单元的通过多点地质统计学随机模拟算法确定初始岩相模型的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度；
步骤608、在当前进行模拟的次数达到所述最大模拟次数时，获取与预先设置的误差容忍度最接近的误差E对应的概率扰动更新后的岩相模型作为所述岩相模型结果，获取与预先设置的误差容忍度最接近的误差E对应的伪储层物性参数模型作为伪储层物性参数模型结果。
此外，所述退火优化单元，具体用于执行：
步骤701、获取待处理的伪储层物性参数模型结果，设置一初始体系动能CΓ(t)和退火终止条件；
步骤702、根据公式计算目标函数E(m(k))；其中，dj为观测到的地震数据；m(k)为第k次迭代的模型参数；
步骤703、根据公式m(k+1)＝m(k)+ξα对所述模型参数进行更新，确定第k+1次迭代的模型参数m(k+1)；ξ为0-1之间的随机数，α为步长；
步骤704、根据公式ΔE＝E(m(k+1))-E(m(k))计算体系势能ΔE；
步骤705、在ΔE<0时，采用m(k+1)替换m(k)；
步骤706、在ΔE≥0时，根据体系温度为T时的概率接受准则ρ(H)＝(e-H/PT)P，在m(k+1)和m(k)中确定概率ρ(H)的较大值对应的模型参数；其中，H＝ΔE+CΓ(t)；P为模型参数总数；
步骤707、在当前迭代的模型参数满足退火终止条件时，将m(k+1)作为储层物性参数反演结果；
步骤708、在当前迭代的模型参数不满足所述退火终止条件时，返回步骤703的模型参数更新过程，直至迭代的模型参数满足所述退火终止条件。
一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现以下步骤：
步骤1、通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数；
步骤2、获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型；
步骤3、统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型；
步骤4、根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型；
步骤5、对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录；
步骤6、将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果；
步骤7、将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
一种计算机设备，包括存储器、处理器及存储在存储上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现以下步骤：
步骤1、通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数；
步骤2、获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型；
步骤3、统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型；
步骤4、根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型；
步骤5、对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录；
步骤6、将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果；
步骤7、将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
本发明实施例提供的一种基于多点地质统计学的地质随机反演方法及装置，将概率扰动反演策略引入到多点地质统计学随机反演中，并利用统计岩石物理和地震正演理论，结合量子退火优化算法形成了一种基于多点地质统计学的地质随机反演方法；本发明能够利用多点地质统计学随机模拟岩相，考虑空间中多点之间的相关性，能够更好的刻画目标体相态；利用概率扰动反演策略，能够减少多点地质统计的模拟次数，减少计算量，在较少的迭代次数中获得随机反演结果，实用性较强。另外通过量子退火优化算法，快速地将产生的具有较差空间连续性的储层物性参数模型优化，从而获得最终的储层物性参数反演结果。本发明能够在较少的迭代次数中同时获得具有较高分辨率和精度的岩相和物性参数的反演结果。随机反演结果能为储层表征和油藏描述提供相对可靠的参考信息。可见，本发明可以解决目前基于多点地质统计学的随机反演方法存在适用范围局限、计算效率低以及仅能获得岩相反演结果的问题。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种基于多点地质统计学的地质随机反演方法的流程图一；
图2为本发明实施例中步骤1的具体实施流程图；
图3为本发明实施例中扫描训练图像获得数据事件的示意图；
图4(a)为本发明实施例中测试模型数据的岩相剖面示意图；
图4(b)为本发明实施例中测试模型数据的孔隙度剖面示意图；
图5为本发明实施例中步骤2的具体实施流程图；
图6(a)为根据概率扰动方法得到的岩相模型的示意图；
图6(b)为根据概率扰动方法得到的岩相模型的孔隙度反演结果的示意图；
图6(c)为根据概率扰动方法得到的岩相模型的孔隙度反演结果误差的示意图；
图7为本发明实施例中步骤3的具体实施流程图；
图8为本发明实施例中的用于多点地质统计学随机反演的一A井数据的示意图；
图9为本发明实施例中的用于多点地质统计学随机反演的一B井数据的示意图；
图10为本发明实施例中步骤5的具体实施流程图；
图11为本发明实施例中的实际地震记录的示意图；
图12(a)为本发明实施例中的实际地震记录的岩相反演结果示意图；
图12(b)为本发明实施例中的实际地震记录的孔隙度反演结果示意图；
图13为本发明实施例中步骤6的具体实施流程图；
图14为本发明实施例中步骤7的具体实施流程图；
图15为本发明实施例中的基于多点地质统计学的地质随机反演装置的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明实施例提供一种基于多点地质统计学的地质随机反演方法，包括：
步骤1、通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数。
步骤2、获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型。
步骤3、统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型。
步骤4、根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型。
步骤5、对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录。
步骤6、将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果。
步骤7、将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
本发明实施例将概率扰动反演策略引入到多点地质统计学随机反演中，并利用统计岩石物理和地震正演理论，结合量子退火优化算法形成了一种基于多点地质统计学的地质随机反演方法；本发明能够利用多点地质统计学随机模拟岩相，考虑空间中多点之间的相关性，能够更好的刻画目标体相态；利用概率扰动反演策略，能够减少多点地质统计的模拟次数，减少计算量，在较少的迭代次数中获得随机反演结果，实用性较强。另外通过量子退火优化算法，快速地将产生的具有较差空间连续性的储层物性参数模型优化，从而获得最终的储层物性参数反演结果。本发明能够在较少的迭代次数中同时获得具有较高分辨率和精度的岩相和物性参数的反演结果。随机反演结果能为储层表征和油藏描述提供相对可靠的参考信息。可见，本发明可以解决目前基于多点地质统计学的随机反演方法存在适用范围局限、计算效率低以及仅能获得岩相反演结果的问题。
为了使本领域的技术人员更好的了解本发明，下面结合上述图1中的步骤1至步骤7对各步骤进行详细阐述：
首先，本发明建立在多点地质统计学模拟的基础上，多点地质统计学以像元为基本模拟单元，利用硬数据作为条件数据，因此模拟结果易于忠实测井数据；同时结合了基于目标的算法的优点，能够较好的再现模拟目标体的结构形态。目前多点地质统计学主要用于离散变量的研究。SNESIM算法是一种适用性较强的多点地质统计学随机模拟方法，计算速度快，效率高，还能有效融合多源信息并能进行多重网格模拟，受到广大地球科学工作者的青睐。SNESIM算法的关键是获取条件概率分布。多点地质统计学中用于表征空间相关性的工具是训练图像，训练图像是一种地质模式的集合，它是能表述实际储层的结构、几何形态及其分布模式的数字化图像，通常被看作是先验的地质概念，主要包含了待模拟区域的主要的特征模式，不要求与精准的测井信息一致或者服从一定的分布。训练图像可以反映地质体空间分布的一般性特征，是特征模式定量化的表示，可以体现出不同特征模式的空间关系。
因此，如图2所示，上述步骤1中通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数，具体可采用如下方式实现：
步骤101、建立研究区域的地质的模拟网格，获得训练图像，将研究区域的地质的条件数据分配至训练图像的网格节点上，并设置一随机路径。
步骤102、根据预先设置的数据模板扫描所述训练图像获得数据事件。
其中所述数据模板为网格节点之间的空间组合关系。
扫描所述训练图像获得数据事件可以如图3所示，其中u为训练图像中的待估计点。另外，如图4(a)和图4(b)所示，其是测试模型数据的岩相(图4(a)所示)和孔隙度(图4(b)所示)剖面，条件数据从中抽取，其分布位置如图中A和B对应的直线所示，该图4(a)和图4(b)中Depth为深度、CDP为深度道集(Common Depth Point)，F为岩相，φ为孔隙度。
步骤103、根据所述数据事件在所述训练图像中检索，得到待估计点处的条件概率分布。
其中，所述待估计点的条件概率分布为：
P{F(u)＝Fk|dn}≈ck(dn)/c(dn)
其中，c(dn)为所述数据事件在训练图像中的有效训练图像中重复的数量；ck(dn)为数据事件的所有重复的数量中F(u)＝Fk的个数；P{F(u)＝Fk|dn}为待估计点u的条件概率分布。
该待估计点的条件概率分布P{F(u)＝Fk|dn}≈ck(dn)/c(dn)是来自如下过程：
在多点地质统计学随机模拟中，待估计点处的条件概率分布函数可以表示为：
其中，该式分母表示数据事件α＝1,…,n出现的概率，相当于数据事件在有效训练图像中重复的数量c(dn)与有效训练图像(如图3中的灰色区域)大小Nn的比值c(dn)/Nn，uα表示数据模板中的其它点；该式的分子表示数据事件和待估计点u取Fk状态同时出现的概率，即数据事件的所有重复个数中F(u)＝Fk的个数ck(dn)和有效训练图像大小Nn的比，可以表示为ck(dn)/Nn，从而待估计点的条件概率分布为：
P{F(u)＝Fk|dn}≈ck(dn)/c(dn)。
步骤104、根据所述待估计点的条件概率分布进行抽样，获得一模拟节点的模拟岩相，作为一次模拟结果，并将所述模拟结果加入到所述条件数据中。
步骤105、随机选取待估计点中的另一模拟节点，重复步骤103至步骤105，直至各待估计点均模拟结束，获得所述研究区域的初始岩相模型。
另外，在本发明实施例中，通过改变随机路径，可以获得其它初始岩相模型。
步骤106、获取预先设置的最大模拟次数。
另外，概率扰动是一种迭代方式的随机反演策略。采用Tau模型，可以将不同来源的约束信息以条件概率的形式融合为一个联合条件概率，运行速度较快、计算量较小。通常模拟几十次左右就可以获得比较满意的反演结果。
因此，如图5所示，上述步骤2中，获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型，具体可以采用如下方式实现：
步骤201、获取预先设置的最大扰动次数和权重因子。
步骤202、对所述初始岩相模型进行概率扰动更新，确定概率扰动更新后的联合概率。
所述概率扰动更新后的联合概率为：P(Fk|dobs)＝(1-r)*i0(u,fk)+r*P(Fk)；其中，为指示变量；P(Fk)表示空间中某未知点u处的岩相分布先验概率；P(Fk|dobs)为该未知点u的概率扰动更新后的联合概率；r表示权重因子，所述权重因子为采用非线性变换的权重因子，可更好地对初始岩相模型进行概率扰动。
步骤203、将所述概率扰动更新后的联合概率P(Fk|dobs)和预先得到的多点概率P(Fk|G)融合，确定扰动融合后的条件概率分布。
所述扰动融合后的条件概率分布为其采用的是Tau模型，其中τ1和τ2为Tau模型参数，a＝(1-P(Fk))/P(Fk)，b＝(1-P(Fk|G))/P(Fk|G)，c＝(1-P(Fk|dobs))/P(Fk|dobs)。
步骤204、将扰动融合后的条件概率分布转化为累积概率分布，根据所述累积概率分布，采用随机抽样方法依次对模拟节点进行抽样，生成概率扰动更新后的岩相模型。
例如根据概率扰动方法可以得到图6(a)的岩相模型，图6(b)的孔隙度反演结果，图6(c)的孔隙度反演结果误差。其中，Δφ为孔隙度反演结果误差。
具体的，如图7所示，上述步骤3中，统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型，可以采用如下方式实现：
步骤301、根据预先获知的测井数据，统计不同岩相中的储层物性参数的分布特征，获得不同岩相中的储层物性参数的分布函数。
步骤302、根据随机抽样方法对概率扰动更新后的岩相模型中的每个网格节点进行随机抽样，获得每个相应网格节点处的物性参数状态值，由每个相应网格节点处的物性参数状态值构成所述伪储层物性参数模型。
由于在转化过程中，没有考虑物性参数变化的空间相关性，仅通过随机抽样获得物性参数分布模型，该模型具有较差的空间相关性，因此称为伪储层物性参数模型。对全部网格节点每执行一次抽样获得一个伪储层物性参数模型，多次随机抽样获得多个伪储层物性参数模型。
另外，在得到伪储层物性参数模型后，根据岩石物理理论即可将其转化为弹性参数(如：速度、密度、波阻抗等)，便于正演。在实际的地震勘探中，地下储层条件十分复杂，不同的孔隙结构、矿物颗粒磨圆度、地层温度、压力条件变化以及泥质含量等因素都会在一定程度上造成弹性参数和储层物性参数岩石物理关系的偏差。确定性的岩石物理模型，没有考虑这种偏差的作用。为了减弱这类影响，在确定性岩石物理关系中加入随机误差项，构成统计岩石物理模型。
因此，上述步骤4中，根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型，可采用如下方式实现：
根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系：e＝fRPM(r)+ε，将所述伪储层物性参数模型转换为弹性参数模型。
其中，fRPM(r)为通过理论岩石物理模型或者经验岩石物理关系获得的确定性岩石物理模型；e为弹性参数；r为储层物性参数；ε为确定性岩石物理模型的误差项，可以通过实际测井资料与确定性岩石物理关系之间的相对差异来估算，常选取均值为零的截断高斯分布，来削弱地下复杂结构对二者关系的影响。
针对不同研究区域的特点，应该选择比较合适的岩石物理模型或统计岩石物理模型以降低在从物性参数到弹性参数转化过程中的不确定性，因此岩石物理模型应谨慎选择。
此外，由于井数据是在深度域获取的(例如图8的A井数据和图9的B井数据，其中φ为孔隙度；Vsh为泥质相对体积；Sw为含水饱和度；Vp为孔隙体积；ρ为密度；Facies为地震相)，前述所产生的岩相模型和伪储层物性参数模型以及弹性参数模型都是以此为条件数据的，因此是深度域的数据，而地震记录是时间域的数据，因此在合成地震数据之间，要将深度域数据转化为时间域数据，在获得时间域的弹性参数模型之后，就可以进行正演获得合成地震记录dsyn，便于与实际地震记录dobs对比，从而选择出能与实际地震记录最匹配的岩相分布和物性参数分布，实现岩相和物性参数的随机反演。
因此，如图10所示，上述步骤5中对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录，可以采用如下方式实现：
步骤501、对所述弹性参数模型进行时深转换获得时间域弹性参数模型。
步骤502、根据所述时间域弹性参数模型计算反射系数。
步骤503、采用褶积模型正演合成地震记录
其中，w为采用褶积模型正演合成的地震记录；b为从实际地震记录中提取的用于合成地震记录的地震子波；R为所述反射系数。
该实际地震记录可以如图11所示，该实际地震记录的岩相反演结果如图12(a)所示，该实际地震记录的孔隙度反演结果如图12(b)所示。
具体的，如图13所示，上述步骤6中，将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果，可以采用如下方式实现：
步骤601、将所述地震记录与预先确定的实际地震记录进行对比，确定所述地震记录与预先确定的实际地震记录的误差E＝‖dobs-dsyn‖2。
其中，dobs为所述实际地震记录；dsyn为正演合成的地震记录。
步骤602、判断所述误差E是否小于或等于预先设置的误差容忍度。
步骤603、在所述误差E小于或等于所述误差容忍度时，将概率扰动更新后的岩相模型作为所述岩相模型结果，将所述伪储层物性参数模型作为伪储层物性参数模型结果。
步骤604、在所述误差E大于所述误差容忍度时，判断进行概率扰动更新的次数是否达到所述最大扰动次数。
步骤605、在进行概率扰动更新的次数未达到所述最大扰动次数时，返回步骤2的进行概率扰动更新的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度。
步骤606、在进行概率扰动更新的次数达到所述最大扰动次数时，判断当前进行模拟的次数是否达到所述最大模拟次数。
步骤607、在当前进行模拟的次数未达到所述最大模拟次数时，修改随机路径，并返回步骤1的通过多点地质统计学随机模拟算法确定初始岩相模型的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度。
步骤608、在当前进行模拟的次数达到所述最大模拟次数时，获取与预先设置的误差容忍度最接近的误差E对应的概率扰动更新后的岩相模型作为所述岩相模型结果，获取与预先设置的误差容忍度最接近的误差E对应的伪储层物性参数模型作为伪储层物性参数模型结果。
在前述阶段所输出的伪储层物性参数模型能够满足事先设定的误差容忍度，使得合成的地震记录与实际地震记录比较接近，与地下实际物性参数分布情况较为接近，但由于其是根据随机抽样获得的，没有考虑空间相关性，因此空间连续性较差，这与实际情况不符，而地震反演通常又具有多解性。因此将该输出的伪储层物性参数模型结果进行优化，提高其空间连续性。量子退火在退火收敛速度和避免陷入局部极小方面有一定优势，对于提高随机反演的精度和效率都有重要价值。
因此，如图14所示，上述步骤7中将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果，可采用如下方式实现：
步骤701、获取待处理的伪储层物性参数模型结果，设置一初始体系动能CΓ(t)和退火终止条件。
步骤702、根据公式计算目标函数E(m(k))；其中，dj为观测到的地震数据；m(k)为第k次迭代的模型参数。
步骤703、根据公式m(k+1)＝m(k)+ξα对所述模型参数进行更新，确定第k+1次迭代的模型参数m(k+1)；ξ为0-1之间的随机数，α为步长。
步骤704、根据公式ΔE＝E(m(k+1))-E(m(k))计算体系势能ΔE。
步骤705、在ΔE<0时，采用m(k+1)替换m(k)。
步骤706、在ΔE≥0时，根据体系温度为T时的概率接受准则ρ(H)＝(e-H/PT)P，在m(k+1)和m(k)中确定概率ρ(H)的较大值对应的模型参数；其中，H＝ΔE+CΓ(t)；P为模型参数总数；C为常数，M为地震道数。
步骤707、在当前迭代的模型参数满足退火终止条件时，将m(k+1)作为储层物性参数反演结果。
步骤708、在当前迭代的模型参数不满足所述退火终止条件时，返回步骤703的模型参数更新过程，直至迭代的模型参数满足所述退火终止条件。
综上所述，本发明实施例提供的基于多点地质统计学的地质随机反演方法具备以下优点：1、应用多点地质统计学随机模拟，考虑空间中多点之间的相关性，更好的表征目标的形态等特征；2、引入概率扰动反演策略，利用Tau模型融合多种来源的约束信息，减少了模拟次数，加快了反演进程，能够提高反演效率，缩短了计算时间；3、本发明结合了量子退火算法，在获得的伪物性参数分布的基础上进行优化，能减小搜索空间，较为快速地获得物性参数的反演结果；4、本发明可以在较少的迭代次数中获得具有较高分辨率的岩相和储层物性参数反演结果，从而为储层评价和油藏的精细描述提供较为可靠的储层模型；5、本发明还可根据需要输出弹性参数的反演结果，便于后续研究；6、在不需要改变反演流程的前提下，利用叠前地震数据和正演算子很容易即可将随机反演从叠后推广到叠前。
另外，如图15所示，本发明实施例还提供一种基于多点地质统计学的地质随机反演装置，包括：
初始岩相模型确定单元11，用于通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数。
概率扰动更新单元12，用于获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型。
伪储层物性参数模型转化单元13，用于统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型。
弹性参数模型转换单元14，用于根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型。
正演合成单元15，用于对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录。
结果确定单元16，用于将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果。
退火优化单元17，用于将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
此外，所述初始岩相模型确定单元11，具体用于执行：
步骤101、建立研究区域的地质的模拟网格，获得训练图像，将研究区域的地质的条件数据分配至训练图像的网格节点上，并设置一随机路径。
步骤102、根据预先设置的数据模板扫描所述训练图像获得数据事件；所述数据模板为网格节点之间的空间组合关系。
步骤103、根据所述数据事件在所述训练图像中检索，得到待估计点处的条件概率分布；所述待估计点的条件概率分布为：
P{F(u)＝Fk|dn}≈ck(dn)/c(dn)
其中，c(dn)为所述数据事件在训练图像中的有效训练图像中重复的数量；ck(dn)为数据事件的所有重复的数量中F(u)＝Fk的个数；P{F(u)＝Fk|dn}为待估计点u的条件概率分布。
步骤104、根据所述待估计点的条件概率分布进行抽样，获得一模拟节点的模拟岩相，作为一次模拟结果，并将所述模拟结果加入到所述条件数据中。
步骤105、随机选取待估计点中的另一模拟节点，重复步骤103至步骤105，直至各待估计点均模拟结束，获得所述研究区域的初始岩相模型。
步骤106、获取预先设置的最大模拟次数。
此外，所述概率扰动更新单元12，具体用于执行：
步骤201、获取预先设置的最大扰动次数和权重因子。
步骤202、对所述初始岩相模型进行概率扰动更新，确定概率扰动更新后的联合概率。
所述概率扰动更新后的联合概率为：P(Fk|dobs)＝(1-r)*i0(u,fk)+r*P(Fk)；其中，为指示变量；P(Fk)表示空间中某未知点u处的岩相分布先验概率；P(Fk|dobs)为该未知点u的概率扰动更新后的联合概率；r表示权重因子，所述权重因子为采用非线性变换的权重因子。
步骤203、将所述概率扰动更新后的联合概率P(Fk|dobs)和预先得到的多点概率P(Fk|G)融合，确定扰动融合后的条件概率分布。
所述扰动融合后的条件概率分布为
其中τ1和τ2为模型参数，a＝(1-P(Fk))/P(Fk)，b＝(1-P(Fk|G))/P(Fk|G)，c＝(1-P(Fk|dobs))/P(Fk|dobs)。
步骤204、将扰动融合后的条件概率分布转化为累积概率分布，根据所述累积概率分布，采用随机抽样方法依次对模拟节点进行抽样，生成概率扰动更新后的岩相模型。
此外，所述伪储层物性参数模型转化单元13，具体用于执行：
步骤301、根据预先获知的测井数据，统计不同岩相中的储层物性参数的分布特征，获得不同岩相中的储层物性参数的分布函数。
步骤302、根据随机抽样方法对概率扰动更新后的岩相模型中的每个网格节点进行随机抽样，获得每个相应网格节点处的物性参数状态值，由每个相应网格节点处的物性参数状态值构成所述伪储层物性参数模型。
此外，所述弹性参数模型转换单元14，具体用于：
根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系：e＝fRPM(r)+ε，将所述伪储层物性参数模型转换为弹性参数模型。
其中，fRPM(r)为通过理论岩石物理模型或者经验岩石物理关系获得的确定性岩石物理模型；e为弹性参数；r为储层物性参数；ε为确定性岩石物理模型的误差项。
此外，所述正演合成单元15，具体用于执行：
步骤501、对所述弹性参数模型进行时深转换获得时间域弹性参数模型。
步骤502、根据所述时间域弹性参数模型计算反射系数。
步骤503、采用褶积模型正演合成地震记录其中，w为采用褶积模型正演合成的地震记录；b为从实际地震记录中提取的用于合成地震记录的地震子波；R为所述反射系数。
此外，所述结果确定单元16，具体用于执行：
步骤601、将所述地震记录与预先确定的实际地震记录进行对比，确定所述地震记录与预先确定的实际地震记录的误差E＝||dobs-dsyn||2；其中，dobs为所述实际地震记录；dsyn为正演合成的地震记录。
步骤602、判断所述误差E是否小于或等于预先设置的误差容忍度。
步骤603、在所述误差E小于或等于所述误差容忍度时，将概率扰动更新后的岩相模型作为所述岩相模型结果，将所述伪储层物性参数模型作为伪储层物性参数模型结果。
步骤604、在所述误差E大于所述误差容忍度时，判断进行概率扰动更新的次数是否达到所述最大扰动次数。
步骤605、在进行概率扰动更新的次数未达到所述最大扰动次数时，返回概率扰动更新单元12的进行概率扰动更新的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度。
步骤606、在进行概率扰动更新的次数达到所述最大扰动次数时，判断当前进行模拟的次数是否达到所述最大模拟次数。
步骤607、在当前进行模拟的次数未达到所述最大模拟次数时，修改随机路径，并返回初始岩相模型确定单元11的通过多点地质统计学随机模拟算法确定初始岩相模型的过程，直至步骤602中的误差E小于或等于预先设置的误差容忍度。
步骤608、在当前进行模拟的次数达到所述最大模拟次数时，获取与预先设置的误差容忍度最接近的误差E对应的概率扰动更新后的岩相模型作为所述岩相模型结果，获取与预先设置的误差容忍度最接近的误差E对应的伪储层物性参数模型作为伪储层物性参数模型结果。
此外，所述退火优化单元17，具体用于执行：
步骤701、获取待处理的伪储层物性参数模型结果，设置一初始体系动能CΓ(t)和退火终止条件。
步骤702、根据公式计算目标函数E(m(k))；其中，dj为观测到的地震数据；m(k)为第k次迭代的模型参数。
步骤703、根据公式m(k+1)＝m(k)+ξα对所述模型参数进行更新，确定第k+1次迭代的模型参数m(k+1)；ξ为0-1之间的随机数，α为步长。
步骤704、根据公式ΔE＝E(m(k+1))-E(m(k))计算体系势能ΔE。
步骤705、在ΔE<0时，采用m(k+1)替换m(k)。
步骤706、在ΔE≥0时，根据体系温度为T时的概率接受准则ρ(H)＝(e-H/PT)P，在m(k+1)和m(k)中确定概率ρ(H)的较大值对应的模型参数；其中，H＝ΔE+CΓ(t)；P为模型参数总数。
步骤707、在当前迭代的模型参数满足退火终止条件时，将m(k+1)作为储层物性参数反演结果。
步骤708、在当前迭代的模型参数不满足所述退火终止条件时，返回步骤703的模型参数更新过程，直至迭代的模型参数满足所述退火终止条件。
值得说明的是，本发明实施例提供的基于多点地质统计学的地质随机反演装置的具体实现方式可以参见上述方法实施例，此处不再赘述。
本发明实施例将概率扰动反演策略引入到多点地质统计学随机反演中，并利用统计岩石物理和地震正演理论，结合量子退火优化算法形成了一种基于多点地质统计学的地质随机反演装置；本发明能够利用多点地质统计学随机模拟岩相，考虑空间中多点之间的相关性，能够更好的刻画目标体相态；利用概率扰动反演策略，能够减少多点地质统计的模拟次数，减少计算量，在较少的迭代次数中获得随机反演结果，实用性较强。另外通过量子退火优化算法，快速地将产生的具有较差空间连续性的储层物性参数模型优化，从而获得最终的储层物性参数反演结果。本发明能够在较少的迭代次数中同时获得具有较高分辨率和精度的岩相和物性参数的反演结果。随机反演结果能为储层表征和油藏描述提供相对可靠的参考信息。可见，本发明可以解决目前基于多点地质统计学的随机反演方法存在适用范围局限、计算效率低以及仅能获得岩相反演结果的问题。
另外，本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现以下步骤：
步骤1、通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数。
步骤2、获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型。
步骤3、统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型。
步骤4、根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型。
步骤5、对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录。
步骤6、将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果。
步骤7、将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
另外，本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现以下步骤：
步骤1、通过多点地质统计学随机模拟算法确定初始岩相模型，并获取预先设置的最大模拟次数。
步骤2、获取预先设置的最大扰动次数和权重因子，并根据概率扰动方法对所述初始岩相模型进行概率扰动更新，生成概率扰动更新后的岩相模型。
步骤3、统计不同岩相中的储层物性参数的分布函数，并根据随机抽样方法将概率扰动更新后的岩相模型转化为伪储层物性参数模型。
步骤4、根据预先设置的统计岩石物理模型中对应研究区域储层物性参数和弹性参数之间的关系，将所述伪储层物性参数模型转换为弹性参数模型。
步骤5、对所述弹性参数模型进行时深转换获得时间域弹性参数模型，确定反射系数，并根据所述反射系数采用褶积模型正演合成地震记录。
步骤6、将所述地震记录与预先确定的实际地震记录进行对比，并根据所述最大扰动次数和最大模拟次数，确定岩相模型结果和伪储层物性参数模型结果。
步骤7、将所述伪储层物性参数模型结果进行量子退火优化，生成储层物性参数反演结果。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。