成矿研究区的成矿模拟方法及装置与流程

1.本发明涉及成矿区确定技术，尤其涉及一种成矿研究区的成矿模拟方法及装置。


背景技术：

2.成矿及成矿作用是一种复杂的动力学过程，矿床成因的基本问题，归根结底是成矿作用的动力学问题。成矿动力学主要是研究成矿作用的速率、机制和过程，是矿床成因的核心问题，传统的矿床学主要总结和归纳已发现矿床的特征、规律和经验，以及对这种规律及事实的解释和进一步的推论，成矿学则是对更大范围的区域内的成矿特征和规律的总结和解释。二者都是立足于经验事实上，其最本质的特征都是以成矿的结果为研究的主要对象，主要目的是发现矿床形成的规律，研究思路则是从结果反推原因，而且这种反推在有很大程度上带有猜想和假想的成分，是非常不严密的。成矿动力学则是以成矿过程为研究的主要对象，目的是发现矿床形成的原因和机制，研究的方法主要是从原因到结果。
3.而近年来随着计算机技术的发展，计算机为成矿动力学提供了很多辅助技术手段，运用计算机技术进行成矿动力学数值模拟成为了可能。综合利用地质、地球物理、地球化学、遥感等研究成果，建立不受时空限制的地质模型，成为现代地球科学研究的重要方法之一。从20世纪九十年代数值模拟技术在地质领域开始广泛应用，地质系统的数值模拟成为当下地质领域的热点和发展趋势，也是地质科学精确化和定量化的必然途径。成矿动力学数值模拟可以模拟和再现成矿的地质过程，并成为现代地球科学研究的重要方法之一。它是在详细的地质工作的基础上，结合地质资料和实验成果，把地质现象定量化，通过归纳、演绎的逻辑推测方法，建立成矿地质过程模型，以帮助人们认识地质构造演化过程。通过数值模拟实验研究，可以深刻理解地质历史时期流体、应力与温度之间的耦合关系；在地质演变中，各历史时期内岩石圈构造、应力应变场、温度场的演化过程，矿体在形成的成矿过程中能量之间的传输与转换及其矿液的聚集过程，矿床的形成过程及其形成机制等。


技术实现要素：

4.有鉴于此，本发明实施例的一个方式是提供一种成矿研究区的成矿模拟方法及装置。
5.本发明一方面提供一种成矿研究区的成矿模拟方法，包括：
6.基于地质的成矿模式与平面、剖面地质资料建立三维地质模型；
7.通过采样实验测试方式至少获得所述地质模型中的数值模拟参数，所述数值模拟参数至少包括以下至少之一：岩石力学参数、构造动力学参数、流体力学参数与热力学参数；
8.将所述三维岩体
‑
构造
‑
通道地质模型导入地质数值模拟参数赋值应用中，进行数值模拟的本构模型的选择与确定，对所述三维岩体
‑
构造
‑
通道地质模型中不同介质材料至少赋予岩石力学参数、构造动力学参数、流体力学参数与热力学参数；
9.根据对成矿研究区实地地质参数统计与样品测试反演，获得成矿过程中最开始的
原始地质参数与最终成矿时各个岩体介质成矿地质参数，基于所述原始地质参数对成矿研究区的初始条件进行赋值；以及根据各个岩体介质成矿地质参数对成矿研究区的边界条件进行设定；
10.以步长设定从小到大的方式，设置数值模拟的步骤，进行动态成矿模拟，并生成相应的模拟动图，当实际的成矿特征与迭代的某个次数相同时，记录其当前参数，在模拟中重新调整迭代次数，使得最终模拟数值模型符合三维成矿模型；
11.数值模拟结束后，导出数值模拟的三维构造应力值、三维流体状态数据；
12.在三维构造应力值、三维流体状态数据中选取符合成矿要求参数，进行布尔运算后即可获得靶区，将所述靶区作为成矿研究区中的成矿区域。
13.优选的，所述基于地质的成矿模式与平面、剖面地质资料建立三维地质模型，包括：
14.基于地质图、巷道工程平面图、钻井编录图、勘探剖面图、地球物理解译剖面图中的至少之一，至少进行矢量化、标准化、要素信息提取的处理，汇总地层、岩性、地球物理性质的地质信息的水平、垂直空间信息，在三维建模应用中进行数据集成，获得地质体的空间位置与接触关系；
15.对地质中的钻孔样、坑道样、地表样的岩石进行近红外光谱分析，获得地质的蚀变类型；进行三维蚀变填图时，基于所述蚀变类型获得蚀变岩性的空间变化信息，指示出热液的主通道位置，作为热液蚀变通道空间位置；
16.对地质中断裂断层进行产状统计，获得成矿期前后的构造运动方向、规模；通过对断层两盘的岩石的地球化学测试，获得其构造运动期次；通过对断层中的矿物进行晶形分析，获得成矿期的构造应力，以此构造断裂空间位置；
17.基于x射线荧光光谱分析、主微量元素的方法获取地质中的携矿岩体、成矿岩体、围岩的钻井岩芯品位数据、坑道探槽品位数据、地表化探数据，以此确定成矿体空间位置；
18.在成矿区板块构造特征、矿集区控矿因素、矿区成矿模式的综合地质约束下，将所述地质体的空间位置与接触关系、所述热液蚀变通道空间位置、所述断裂空间位置及所述矿体空间位置进行综合，在三维建模应用的约束下，建立三维岩体
‑
构造
‑
通道地质模型。
19.优选的，所述岩石力学参数包括以下至少之一：密度、剪切模量、体积模量、粘聚力、内摩擦角、膨胀角和抗拉强度。
20.优选的，所述构造动力学参数包括以下至少之一：岩体运动速度、构造应力；
21.对应地，所述获得所述地质模型中的数值模拟参数，包括：
22.通过对成矿区的断裂断层产状统计获得成矿期前后的构造运动方向、规模距离；通过对断层两盘的岩石的地球化学测试，获得成矿研究区构造运动期次及时间；以及，通过对断层中的矿物晶形分析，获得成矿研究区成矿期的构造应力；
23.综合构造运动方向、规模距离，以及与运动时间获得成矿研究区运动平均速度。
24.优选的，所述方法还包括：
25.对岩石包裹体进行测试，获得成矿期流体的流体压力与流体温度，对围岩岩石样品进行温度渗透实验，获得岩石的渗透率与孔隙率，以及岩石热传导数据。
26.本发明的另一方面提供一种成矿研究区的成矿模拟装置，包括：
27.建立单元，用于基于地质的成矿模式与平面、剖面地质资料建立三维地质模型；
28.获得单元，用于通过采样实验测试方式至少获得所述地质模型中的数值模拟参数，所述数值模拟参数至少包括以下至少之一：岩石力学参数、构造动力学参数、流体力学参数与热力学参数；
29.赋值单元，用于将所述三维岩体
‑
构造
‑
通道地质模型导入地质数值模拟参数赋值应用中，进行数值模拟的本构模型的选择与确定，对所述三维岩体
‑
构造
‑
通道地质模型中不同介质材料至少赋予岩石力学参数、构造动力学参数、流体力学参数与热力学参数；
30.设定单元，用于根据对成矿研究区实地地质参数统计与样品测试反演，获得成矿过程中最开始的原始地质参数与最终成矿时各个岩体介质成矿地质参数，基于所述原始地质参数对成矿研究区的初始条件进行赋值；以及根据各个岩体介质成矿地质参数对成矿研究区的边界条件进行设定；
31.模拟单元，用于以步长设定从小到大的方式，设置数值模拟的步骤，进行动态成矿模拟，并生成相应的模拟动图，当实际的成矿特征与迭代的某个次数相同时，记录其当前参数，在模拟中重新调整迭代次数，使得最终模拟数值模型符合三维成矿模型；
32.导出单元，用于在数值模拟结束后，导出数值模拟的三维构造应力值、三维流体状态数据；
33.成矿确定单元，用于在三维构造应力值、三维流体状态数据中选取符合成矿要求参数，进行布尔运算后即可获得靶区，将所述靶区作为成矿研究区中的成矿区域。
34.优选的，所述建立单元，还用于：
35.基于地质图、巷道工程平面图、钻井编录图、勘探剖面图、地球物理解译剖面图中的至少之一，至少进行矢量化、标准化、要素信息提取的处理，汇总地层、岩性、地球物理性质的地质信息的水平、垂直空间信息，在三维建模应用中进行数据集成，获得地质体的空间位置与接触关系；
36.对地质中的钻孔样、坑道样、地表样的岩石进行近红外光谱分析，获得地质的蚀变类型；进行三维蚀变填图时，基于所述蚀变类型获得蚀变岩性的空间变化信息，指示出热液的主通道位置，作为热液蚀变通道空间位置；
37.对地质中断裂断层进行产状统计，获得成矿期前后的构造运动方向、规模；通过对断层两盘的岩石的地球化学测试，获得其构造运动期次；通过对断层中的矿物进行晶形分析，获得成矿期的构造应力，以此构造断裂空间位置；
38.基于x射线荧光光谱分析、主微量元素的装置获取地质中的携矿岩体、成矿岩体、围岩的钻井岩芯品位数据、坑道探槽品位数据、地表化探数据，以此确定成矿体空间位置；
39.在成矿区板块构造特征、矿集区控矿因素、矿区成矿模式的综合地质约束下，将所述地质体的空间位置与接触关系、所述热液蚀变通道空间位置、所述断裂空间位置及所述矿体空间位置进行综合，在三维建模应用的约束下，建立三维岩体
‑
构造
‑
通道地质模型。
40.优选的，所述岩石力学参数包括以下至少之一：密度、剪切模量、体积模量、粘聚力、内摩擦角、膨胀角和抗拉强度。
41.优选的，所述构造动力学参数包括以下至少之一：岩体运动速度、构造应力；
42.对应地，所述获得单元，还用于：
43.通过对成矿区的断裂断层产状统计获得成矿期前后的构造运动方向、规模距离；通过对断层两盘的岩石的地球化学测试，获得成矿研究区构造运动期次及时间；以及，通过
对断层中的矿物晶形分析，获得成矿研究区成矿期的构造应力；
44.综合构造运动方向、规模距离，以及与运动时间获得成矿研究区运动平均速度。
45.优选的，所述获得单元，还用于：
46.对岩石包裹体进行测试，获得成矿期流体的流体压力与流体温度，对围岩岩石样品进行温度渗透实验，获得岩石的渗透率与孔隙率，以及岩石热传导数据。
47.本发明的另一方面提供一种电子设备，包括：处理器和用于存储处理器可执行指令的存储器，其中，所述处理器被配置为在调用存储器中的可执行指令时，能够执行所述的成矿研究区的成矿模拟方法的步骤。
48.本发明的另一方面提供一种计算机可读存储介质，其上存取有计算机指令，其特征在于，所述指令被处理器执行时实现所述的成矿研究区的成矿模拟方法的步骤。
49.本发明实施例的成矿研究区的成矿模拟方法及装置，通过对成矿研究区域的地层、岩体、构造断裂等进行精细三维建模，获得成矿关键期的动力学构造流体通道模型，借助于传统矿床研究方法，获得成矿流体的温度、流体压力及动力学的运动方向，应力大小等数据状态，在经典成矿理论的指导下，对成矿区域、深部剖面、矿床三方面进行成矿动力学数值模拟，揭示矿体与各地质要素的空间分布规律及控矿特征分析，了解成矿机制与演化过程，进一步验证及细化成矿理论。本发明是将复杂的成矿模拟模式流程化规范化，将复杂的多个过程转化为流程化的建模
‑
模拟过程，在前期准备好相应的数据后，可以按照要求填入加入模拟调整参数，迭代调整即可获得模拟结果，该过程大大降低了建模与模拟的复杂性，快捷地实现了空间
‑
时间维度的四维成矿模拟成矿与预测的结果。
附图说明
50.图1为本发明实施例的成矿研究区的成矿模拟方法的流程示意图；
51.图2为本发明实施例的成矿相关要素的三维建模过程示意图；
52.图3为本发明实施例的地质模型中的数值模拟参数采集示意图；
53.图4为本发明实施例的成矿研究区的成矿模拟装置的组成结构示意图。
具体实施方式
54.以下结合附图，详细阐明本发明技术方案的实质。
55.图1为本发明的成矿研究区的成矿模拟方法的流程示意图，如图1所示，本发明实施例的成矿研究区的成矿模拟方法包括以下处理步骤：
56.步骤101，基于地质的成矿模式与平面、剖面地质资料建立三维地质模型。
57.本发明实施例中的建立三维地质模型，具体包括：
58.基于地质图、巷道工程平面图、钻井编录图、勘探剖面图、地球物理解译剖面图中的至少之一，至少进行矢量化、标准化、要素信息提取的处理，汇总地层、岩性、地球物理性质的地质信息的水平、垂直空间信息，在三维建模应用中进行数据集成，获得地质体的空间位置与接触关系；
59.对地质中的钻孔样、坑道样、地表样的岩石进行近红外光谱分析，获得地质的蚀变类型；进行三维蚀变填图时，基于所述蚀变类型获得蚀变岩性的空间变化信息，指示出热液的主通道位置，作为热液蚀变通道空间位置；
60.对地质中断裂断层进行产状统计，获得成矿期前后的构造运动方向、规模；通过对断层两盘的岩石的地球化学测试，获得其构造运动期次；通过对断层中的矿物进行晶形分析，获得成矿期的构造应力，以此构造断裂空间位置；
61.基于x射线荧光光谱分析、主微量元素的方法获取地质中的携矿岩体、成矿岩体、围岩的钻井岩芯品位数据、坑道探槽品位数据、地表化探数据，以此确定成矿体空间位置；
62.在成矿区板块构造特征、矿集区控矿因素、矿区成矿模式的综合地质约束下，将所述地质体的空间位置与接触关系、所述热液蚀变通道空间位置、所述断裂空间位置及所述矿体空间位置进行综合，在三维建模应用的约束下，建立三维岩体
‑
构造
‑
通道地质模型。
63.步骤102，通过采样实验测试方式至少获得所述地质模型中的数值模拟参数。
64.其中，所述数值模拟参数至少包括以下至少之一：岩石力学参数、构造动力学参数、流体力学参数与热力学参数。所述岩石力学参数包括以下至少之一：密度、剪切模量、体积模量、粘聚力、内摩擦角、膨胀角和抗拉强度。
65.所述构造动力学参数包括以下至少之一：岩体运动速度、构造应力；对应地，所述获得所述地质模型中的数值模拟参数，包括：通过对成矿区的断裂断层产状统计获得成矿期前后的构造运动方向、规模距离；通过对断层两盘的岩石的地球化学测试，获得成矿研究区构造运动期次及时间；以及，通过对断层中的矿物晶形分析，获得成矿研究区成矿期的构造应力；综合构造运动方向、规模距离，以及与运动时间获得成矿研究区运动平均速度。
66.本发明实施例的获取数值模拟参数还包括：对岩石包裹体进行测试，获得成矿期流体的流体压力与流体温度，对围岩岩石样品进行温度渗透实验，获得岩石的渗透率与孔隙率，以及岩石热传导数据。
67.步骤103，根据对成矿研究区实地地质参数统计与样品测试反演，获得成矿过程中最开始的原始地质参数与最终成矿时各个岩体介质成矿地质参数，基于所述原始地质参数对成矿研究区的初始条件进行赋值；以及根据各个岩体介质成矿地质参数对成矿研究区的边界条件进行设定。
68.步骤104，以步长设定从小到大的方式，设置数值模拟的步骤，进行动态成矿模拟，并生成相应的模拟动图，当实际的成矿特征与迭代的某个次数相同时，记录其当前参数，在模拟中重新调整迭代次数，使得最终模拟数值模型符合三维成矿模型。
69.步骤105，数值模拟结束后，导出数值模拟的三维构造应力值、三维流体状态数据。
70.步骤106，在三维构造应力值、三维流体状态数据中选取符合成矿要求参数，进行布尔运算后即可获得靶区，将所述靶区作为成矿研究区中的成矿区域。
71.以下通过具体示例，进一步阐明本发明实施例的技术方案的实质。
72.本发明实施例的目的就是通过对成矿研究区域的地层、岩体、构造断裂等进行精细三维建模，获得成矿关键期的的动力学构造流体通道模型，并借助于传统矿床研究，获得成矿流体的温度、流体压力及动力学的运动方向，应力大小等数据状态，在经典成矿理论的指导下，对成矿区域、深部剖面、矿床三方面进行成矿动力学数值模拟，揭示矿体与各地质要素的空间分布规律及控矿特征分析，了解成矿机制与演化过程，进一步验证及细化成矿理论。对成矿模拟的结果进行后处理，也能为深部的隐伏矿体预测提供重要参考依据。
73.本发明实施例需要确定以下预设条件：
74.矿床成矿类型为热液成矿，有明确的成矿作用模式，在成矿期有明确的动力学构
造运动或者热液流体形成的热固耦合成矿作用；可以根据地质资料建立成矿关键期的地层
‑
岩体空间特征，对断裂构造的空间形态有详细的资料，通过平面与剖面地质数据，如探槽露头、地球物理资料，钻井及勘探线数据、井下采矿剖面等建立精准的三维地质模型；通过地球化学研究获得成矿期的主要动力学参数如构造运动的方向，应力大小，和成矿热液的温度、流体压力、流速以及围岩的岩石各种动力学物理性质，初始状态与临界状态等。
75.具体地，本发明实施例的成矿研究区的成矿模拟方法包括以下处理步骤：
76.通过成矿模式与平面、剖面地质资料建立三维地质模型；
77.通过测试成矿研究区的成矿期包裹体样品，获得成矿期的温度、压力等成矿临界状态；
78.通过测试围岩与赋矿岩体的岩石力学特性，获得数值模拟的基本参数；
79.根据成矿模式与成矿区研究，对三维地质模型赋予初始状态条件；
80.调整参数，借助临界条件数据进行模拟，直到获得与临界条件相当的模拟模型；
81.成矿过程在模拟中可以获得；将临界模型进行三维后处理，获得三维成矿范围与成矿潜力区。
82.图2为本发明的成矿相关要素的三维建模过程示意图，如图2所示，本发明实施例的三维建模过程是基于两种主流建模方法，即数据驱动建模与知识框架建模，将成矿研究区的地质信息数据转化为数据模型，在矿床成矿模式的地质约束下，进行三维地质模型的构建，具体步骤如下：
83.首先是三维地质框架信息的收集。如二维的地质图，巷道工程平面图、钻井编录图、勘探剖面图、地球物理(重电磁震)解译剖面图，经过autocad、mapgis、arcgis等多种软件的矢量化、标准化、要素信息提取，将成矿研究区的地层、岩性、地球物理性质等地质信息的水平、垂直空间信息进行汇总，在三维建模软件(如gocad)中进行数据集成，获得地质体的空间位置与接触关系；
84.其次是携矿岩体、成矿岩体、围岩的空间位置确定。主要方法是借助x射线荧光光谱分析(xrf)、主微量元素分析等方法测试成矿元素的品位。测试的主要对象为钻井岩芯品位数据、坑道探槽品位数据、地表化探数据等。通过取样点的空间位置信息与数值，借助地质统计学与三维建模软件插值功能，可以建立已知矿体的三维模型。如果矿山已有矿体三维模型，则可以借鉴使用；
85.热液通道主要是通过不同的岩性蚀变进行通道位置的指示。岩性蚀变在地下坑道中可以通过穿脉目视解译进行追溯，但效果最好的是近红外光谱蚀变解译。通过对钻孔样、坑道样、地表样的岩石近红外光谱分析，获得蚀变类型。在三维蚀变填图的过程中获得蚀变岩性的空间变化信息，进一步指示了热液的主通道位置；
86.构造运动则是通过对成矿区的断裂断层产状统计获得成矿期前后的构造运动方向、规模；通过对断层两盘的岩石的地球化学测试，如锆石测年，获得其构造运动期次；通过对断层中的矿物晶形分析，获得其成矿期的构造应力；
87.最终在成矿区板块构造特征、矿集区控矿因素、矿区成矿模式的综合地质约束下，将以上地质信息进行综合，在三维建模软件(gocad)的约束下，建立三维岩体
‑
构造
‑
通道地质模型。
88.图3为本发明实施例的地质模型中的数值模拟参数采集示意图，如图3所示，本发
明实施例的值模拟参数主要包括岩石力学参数，构造动力学参数，流体力学参数与热力学参数。以上参数均通过采样实验测试获得参考值。
89.岩石力学参数包括密度、剪切模量、体积模量、粘聚力、内摩擦角、膨胀角与抗拉强度。以上测试分别可以根据三维模型中的不同岩性分别取样，在三轴应力测试的仪器下进行多次测试，获取平均的参考值。
90.构造动力学参数包括岩体运动速度、构造应力。岩体的成矿期构造运动可以通过对成矿区的断裂断层产状统计获得成矿期前后的构造运动方向、规模距离；通过对断层两盘的岩石的地球化学测试，如锆石测年，获得其构造运动期次及时间；通过对断层中的矿物晶形分析，获得其成矿期的构造应力；综合相对运动距离获得与运动时间获得其运动平均速度。
91.本发明实施例中，构造运动速度
‑
构造地质信息统计计算如下：
92.利用构造地质学知识，在研究区进行构造信息踏勘统计，将矿床内破碎带构造变形具有与热液蚀变相似的分带特征，断裂带主断裂面和次级断裂面及其线理、(擦痕)产状、规模、运动方式等信息做出统计后，计算构造玫瑰图，反演获得断层的构造运动方向、运动速度等信息。
93.基于前述信息构造应力
‑
石英变形和电子背散射衍射(ebsd)组构。具体地，采用石英动态重结晶亚颗粒大小法对新城金矿床蚀变带差应力值进行计算，金属和矿物的实验研究均表明，动态重结晶颗粒的大小与形变达到稳态时的应力大小有关：应力越大，颗粒越细。采用线截法，在autocad中对新城金矿床不同构造变形期样品动态重结晶颗粒的大小完成统计。利用关系式：σ1‑
σ3＝ad
‑
m
进行测算。其中，d为动态重结晶颗粒大小(μm)，a为6.1，m为0.68。一般动态重结晶颗粒的大小在几微米至几十微米范围。关系式计算得到新城金矿床成矿前差应力值变化范围为85.63～106.85mpa，成矿期差应力变化范围为65.91～76.77mpa。
94.对岩石包裹体做测试，获得成矿期流体的流体压力与流体温度，对围岩岩石样品进行温度渗透实验获得岩石的渗透率与孔隙率与岩石热传导数据。
95.利用成矿岩体中的石英包裹体，通过对其薄片的冷热处理，观察固液气相对象形态变化，反演出其形成时的温度。
96.本发明实施例中，设包裹体的前提为：均匀体系。包裹体形成时，被捕获的流体是均匀体系，即主矿物是在均匀体系中生长的。封闭体系。充填(滞留)在晶体缺陷中的流体为主矿物封闭，形成独立的封闭体系，没有外来物质的加入和内部物质的逸出。等容体系。包裹体形成后，体积基本恒定不变，保持等容体系的特点，因而可以利用各种与之有关的物理化学相图。
97.可以通过以下方式确定石英包裹体：
98.1、冷冻法：指在包裹体冷却到室温以下时观察液相向固相转变(即固化)过程。基本原理是通过在冷台上改变温度，观察包裹体所发生的相变过程。
99.符合拉乌尔定律—对于稀浓度溶液而言，溶液的冰点下降数值与溶质的种类及性质无关，而仅仅取决于溶解在水(溶剂)中的溶质的浓度；对于具有相同浓度的各种溶质，其冰点的下降温度也相同。
100.2、均一法：根据包裹体的基本假设和前提，包裹体所捕获的流体为原始均匀的单
一相流体，它们充满着整个包裹体空间。随着温度下降，流体(气体或液体)的收缩系数大于固体(主矿物)的收缩系数，包裹体将沿着等容线演化，一直到两相界面的位置，如果原来捕获的是大于临界密度的流体，则分离出一个气相，气体逸出后，由于表面张力的影响，气体在有利位置形成球形的气泡；如果原来捕获的是小于临界密度的富气体流体，则气体在流体中凝聚出一个液相，形成具有一个大气泡的两相包裹体。
101.如果在冷热台中升温，则可看到可逆的相变化现象：首先看到的是随着温度的升高气、液相的比例发生变化，而当升到一定温度时，就发生了相的转变，即从两相(或多相)转变成一个相，也即达到了相的均一，这时的温度即为均一温度。
102.本发明实施例中，测温分析的原理实现为，只要在光学显微镜上附加一种测温设备，就能在地质上有意义的各种透明(或半透明)矿物中广泛应用。该方法是在详细观察和辨认包裹体中含流体的各种物相(固相、气相、液相)基础上，通过升温或冷冻来测量各种瞬间相变化的温度。相较于适用0～1500k测温范围的linkam
‑
ts150型热台，linkam thmsg
‑
600型冷热台可以实现测温分析，从而确定出包裹体。
103.以下说明渗透率的确定方式。
104.根据达西定律q/s＝
‑
k
△
p/ηl，式中，q为流量(m3/s)；s为样品横截面积(m2)；l为样品长度(m)；η为流体黏滞系数(pa
·
s)；k为渗透率(m2)；δp为样品上、下游的压力差(pa)。在岩样的上、下游端施加稳定的压力差δp，通过测量流经样品的流量q得到渗透率，或者保持恒定的流量q而测量上、下游端的压力差δp而得到渗透率。
105.本发明实施例中，对超低渗透率测量还可以采用周期振荡法，其测量下限可达10
‑
10
μm 2。周期振荡法最早由kranz等借鉴热扩散系数测量而提出，被运用到测量低渗岩石渗透率上。具体地，将测试岩芯作为一个衰减与阻尼器。首先在岩芯的一端输入一个振幅及频率一致的正弦振荡压力波，穿过岩芯在另一端输出一个振幅及相位发生衰减和延迟的同频率正弦波，振幅衰减和相位延迟与岩石渗透率有关。正弦压力波在岩芯中的传播过程类似一维扩散：
[0106][0107]
式(1)中：p为孔隙压力；k为渗透率；μ为流体的黏滞系数；β
s
为岩芯
‑
流体系统的比储流率。式(1)的初始条件和边界条件为：
[0108][0109]
式(2)中：x＝0选择在下游端面处；sd为下游容器体积结合初始条件和边界条件式(2)式(1)的解为
[0110]
[0111]
从式(3)可以看出，周期振荡法可看成稳态法和脉冲衰减法相结合的一种混和法。下游压力对上游应力的响应由两部分组成，既存在由于压力突然
[0112]
变为正弦波的一个瞬时响应(也呈指数衰减)，也存在一个稳定的正弦周期引起的稳态响应从方程(3)的前半部分看出，相对上游的压力，下游的响
[0113]
应会在幅值上发生衰减(衰减因子α)和相位发生延迟θ。通过测量上、下游幅值比α和相位延迟θ即可计算样品的渗透率。
[0114]
根据fischer研究结论，上、下游压力幅值比α和相位延迟θ可以表示为2个无量纲参数ψ和γ的函数：
[0115][0116]
式(4)中：
[0117][0118]
式(5)中：γ、ψ为叠加过程中产生的2个无量纲值。由以上的各类参数得到的渗透率为
[0119][0120]
式(6)中：t为孔隙压力的振荡周期(s)。
[0121]
由以上式子可知，从实验室中测得的振幅比α和相位移θ，求出了2个无量纲值ψ和γ，并由式(6)可得渗透率k。
[0122]
本发明实施例中，岩石传热导数(热导率)通过以下方式确定。
[0123]
常温常压下测量热导率有很多办法，包括稳态法和瞬态法两个大类，热流计法、保护热板法、圆管法、热线法、闪光法、瞬变平面热源法。理想的热传导系数测量方法应该在原位进行，但是由于实验室测量方法简便，测量精度高，所以大量的岩石热传导系数测量是在实验室进行。在实验室测定时，应该使岩石样品处于原始湿度与温度的状态下进行。
[0124]
稳态法指的是实验中待测试样上温度分布达到稳定后进行测量，其分析的出发点是稳态的导热微分方程，能直接测得导热系数。其特点是实验公式简单，需要测量导热量(直接或间接地)和若干点的温度。缺点就是测试周期太长，装样复杂，比较难以应用于需要测试大量样品的实验中，因此使用该类方法的人并不多。
[0125]
与稳态法对应的是瞬态法，也是比较常见的热导率测量方法。指测量过程中试样温度随时间变化，需要测量试样上若干点的温度随时间的化。其分析的出发点是瞬态导热微分方程，特点是公式复杂，但是样品组装简单，试验周期短.瞬态法可以同时测定样品的
热导率和热扩散系数而且对样品温度的边界条件要求不高。
[0126]
测量原理是对处于热平衡状态的试样施加某种热干扰，同时测量试样对热干扰的响应(温度或热流随时间的变化)，然后根据响应曲线确定热物性参数的数值。
[0127]
运用热常数分析仪测量热传导系数，需要假设探头位于无限大样品内，推导出其热导方程并求解。
[0128]
当热常数分析仪探头被电加热时，可得电阻升高随时间的方程如下(1
‑
1)：
[0129]
r＝r0{1 α[δt
i

△
t(τ)]}
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1
‑
1)
[0130]
式中：r——传感器电阻，r0是传感器被加热前的电阻；
[0131]
α——电阻温度系数(tcr)；
[0132]
δt
i
——镍与绝缘薄层的温度差分；
[0133]
δt(τ)——样品表面温度升高值。
[0134]
从方程(2
‑
1)得到传感器记录的温度增加值如公式(1
‑
2)所示：
[0135][0136]
δt
i
实际上表示了传感器记录温升值与样品表面温升值的差值。
[0137]
温度差分δt
i
在极短时间δt
i
后变为常数，可以估计如(1
‑
3)所示如下：
[0138][0139]
式中：δ——绝缘层厚度，mm；
[0140]
k
i
——绝缘层材料的热扩散系数，mm2/s。
[0141]
样品表面温度升高与时间函数的关系式如(1
‑
4)所示：
[0142][0143]
式中：p0——传感器加热功率，w；
[0144]
a——探头的半径，mm；
[0145]
λ——测试样品的热传导系数，w/(m.k)；
[0146]
d(τ)——与尺寸无关的时间依赖方程，如(1
‑
5)。
[0147][0148]
式中：t——瞬态记录对应的时间，s；
[0149]
θ——特征时间，定义为(1
‑
6)。
[0150][0151]
由式(1
‑
4)可以看出试样表面的平均温升δt(τ)与d(τ)呈线性关系，其中为直线的斜率，通过斜率即可求得试样的导热系数。通过多次试探k值，使得平均温升δt(τ)与d(τ)呈严格线性对应关系，就可以得到热传导系数。再通过导热系数、热扩散系数、体
积比热之间的关系求出体积比热。
[0152]
在对三维模型进行数值模拟前，对岩石等介质进行属性赋值，以某金矿床为例，根据以上方法测试，获得的岩石力学参数如下表1所示：
[0153][0154]
表1
[0155]
根据实地统计与样品测试反演，获得成矿过程中最开始的原始地质参数与最终成矿时各个岩体介质成矿地质参数。
[0156]
上述金矿床的结构分析和显微组织研究表明，在成矿前，北西
‑
南东(nw
‑
se)在挤压条件下经历了韧性
‑
脆性左线剪切变形，差异应力为61.37～111.09mpa。应用nw
‑
se的胶家金矿带北段前述金矿床的有限元数值模拟，主压应力为100mpa，ne
‑
sw剪切力为40mpa。如表2所示。
[0157][0158]
表2
[0159]
上述金矿床的成矿石英包裹体压力范围为1082
‑
3444，1041
‑
2609bar，848
‑
2022bar，1852
‑
3375。取其最大值作为模拟的新城的最小孔隙压力值，即337500000，8.4e7
‑‑‑‑
3.375e8pa。
[0160]
以上测试数据即可作为数值模拟的参数、初始条件与边界条件参与到运算中。
[0161]
本发明实施例中，将三维地质建模软件已经建立的块体模型导入到数值模拟软件中，以flac3d软件为例，导入inp格式的块体文件，其代码如下：
[0162]
model new
[0163]
model configure thermal fluid
[0164]
zone import'gjl.inp'
[0165]
zone group'gjl'
[0166]
zone import'll.inp'
[0167]
zone group'll'range group'gjl'not
[0168]
zone import'fault.inp'
[0169]
zone group'fault'range group'll'not group'gjl'not
[0170]
zone import'fault
‑
b.inp'
[0171]
zone group'fault
‑
b'range group'fault'not group'll'not group'gjl'not
[0172]
zone import'fault
‑
t.inp'
[0173]
zone group'fault
‑
t'range group'fault
‑
b'not group'fault'not group'll'not group'gjl'not
[0174]
zone attach by
‑
face。
[0175]
本发明实施例中，对于成矿数值模拟，以岩石作为典型的应用对象，所以选择典型的摩尔库伦模型(zone cmodel assign mohr
‑
coulomb)。
[0176]
本发明实施例中，对模型中不同介质材料赋予力学参数、流体力学参数等，根据已获得属性值对不同介质赋予相应的参数。示例如下：
[0177]
zone initialize density 2680range group'll'
[0178]
zone initialize density 2650range group'gjl'
[0179]
zone initialize density 2810range group'fault
‑
b'
[0180]
zone initialize density 2810range group'fault
‑
t'
[0181]
zone initialize density 2000range group'fault'
[0182]
zone prop bulk 52e9 shear 34e9 coh 5.3e7 fri 52ten 7.5e6 dil 2range group'll'
[0183]
zone prop bulk 55e9 shear 38e9 coh 5.4e7 fri 53ten 8.2e6 dil 2range group'gjl'
[0184]
；zone prop bulk 52e9 shear 34e9 coh 5.3e7 fri 52ten 7.5e6 dil 2range group'fault
‑
t'
[0185]
；zone prop bulk 55e9 shear 38e9 coh 5.4e7 fri 53ten 8.2e6 dil 2range group 'fault
‑
b'
[0186]
zone prop bulk 57e9 shear 41e9 coh 4.2e7 fri 54ten 7.6e6 dil 2range group'fault
‑
b'
[0187]
zone prop bulk 57e9 shear 41e9 coh 4.2e7 fri 54ten 7.6e6 dil 2range group 'fault
‑
t'
[0188]
zone prop bulk 9e9 shear 4e9 coh 0.7e7 fri 40ten 0.9e6 dil 2range group 'fault'
[0189]
zone fluid cmodel assign isotropic
[0190]
zone fluid property permeability 5e
‑
12range group'll'
[0191]
zone fluid property porosity 0.5range group'll'
[0192]
zone fluid property permeability 5e
‑
12range group'gjl'
[0193]
zone fluid property porosity 0.5range group'gjl'。
[0194]
本发明实施例中，还对初始条件进行赋值，如对构造应力、运动速度、空隙压力、流体温度进行赋值等，根据获得的初始条件，对模型的对应对象进行赋值，如构造应力赋值、温度场赋值、流体渗流初始压力赋值。示例如下：
[0195]
zone fluid cmodel assign anisotropic
[0196]
zone fluid property permeability
‑
xz 1e
‑
6range group'fault
‑
t'
[0197]
zone fluid property porosity 0.9range group'fault
‑
t'
[0198]
zone fluid property permeability
‑
xz 1e
‑
2range group'fault
‑
b'
1000range group'fault
‑
b'
[0229]
zone thermal property conductivity 1.1expansion 3e
‑
6specific
‑
heat 1000range group'fault
‑
t'
[0230]
zone thermal property conductivity 0.7expansion 8e
‑
6specific
‑
heat 1000range group'fault'；
[0231]
zone gridpoint initialize temperature 20
[0232]
zone face apply temperature 5000range group'fault
‑
b'position
‑
z
‑
1501
‑
1499。
[0233]
赋值完成后，开始模拟。具体地，设定好数值模拟的步骤，其步长设定从小到大，实现动态成矿模拟动图，当实际的成矿特征与迭代的某个次数相同时，记录其当前参数，在模拟中重新调整迭代次数，使得最终模拟数值模型符合三维成矿模型。记录结果示例如下：
[0234]
model fluid active on
[0235]
model thermal active on
[0236]
model mechanical active on；
[0237]
model thermal timestep fix 6.48e3
[0238]
model gravity 10
[0239]
model solve ratio 1e
‑5[0240]
；model solve cycles 5000。
[0241]
本发明实施例中，模拟结束后，导出模拟的三维构造应力值、三维流体状态数据进行后处理。将数值模拟结果利用list命令进行导出：zone gridpoint list position导出各个数据的位置信息(xyz坐标)，zone list stress导出压力数据。流体饱和度，流体压力等参数同样可以利用list导出，导出数据为txt格式。将该数据调整为gocad数据导入，选取符合成矿要求参数，布尔运算后即可获得靶区。
[0242]
本发明实施例中，通过将导出的模拟数据导入三维模型，与现在矿体三维模型作对比，不符合则调整模拟迭代次数，直至模型符合为止；当模拟模型与实际矿体模型相符时，根据成矿控矿因素，如构造破碎带(构造变性模拟模型)、流体渗透压模型等输入到三维地质模型中进行后处理，在合适的阈值约束下，获得有利的成矿预测区，进行布尔运算即可获得最佳成矿预测区。
[0243]
图4为本发明实施例的成矿研究区的成矿模拟装置的组成结构示意图，如图4所示，本发明实施例的成矿研究区的成矿模拟装置包括：
[0244]
建立单元40，用于基于地质的成矿模式与平面、剖面地质资料建立三维地质模型；
[0245]
获得单元41，用于通过采样实验测试方式至少获得所述地质模型中的数值模拟参数，所述数值模拟参数至少包括以下至少之一：岩石力学参数、构造动力学参数、流体力学参数与热力学参数；
[0246]
赋值单元42，用于将所述三维岩体
‑
构造
‑
通道地质模型导入地质数值模拟参数赋值应用中，进行数值模拟的本构模型的选择与确定，对所述三维岩体
‑
构造
‑
通道地质模型中不同介质材料至少赋予岩石力学参数、构造动力学参数、流体力学参数与热力学参数；
[0247]
设定单元43，用于根据对成矿研究区实地地质参数统计与样品测试反演，获得成矿过程中最开始的原始地质参数与最终成矿时各个岩体介质成矿地质参数，基于所述原始
processor)、应用专用集成电路(asic，application specific integrated circuit)、dsp、可编程逻辑器件(pld，programmable logic device)、复杂可编程逻辑器件(cpld，complex programmable logic device)、现场可编程门阵列(fpga，field
‑
programmable gate array)、通用处理器、控制器、微控制器(mcu，micro controller unit)、微处理器(microprocessor)、或其他电子元件实现，也可以结合一个或多个射频(rf，radio frequency)天线实现，用于执行前述实施例的网络数据收集方法的步骤。
[0265]
在本公开实施例中，图4示出的成矿研究区的成矿模拟装置中各个单元执行操作的具体方式已经在有关该方法的实施例中进行了详细描述，此处将不做详细阐述说明。
[0266]
本发明实施例还记载了提供一种电子设备，包括：处理器和用于存储处理器可执行指令的存储器，其中，所述处理器被配置为在调用存储器中的可执行指令时，能够执行前述实施例的成矿研究区的成矿模拟方法的步骤。
[0267]
本发明实施例还记载了一种计算机可读存储介质，其上存取有计算机指令，其特征在于，所述指令被处理器执行时实现前述实施例的成矿研究区的成矿模拟方法的步骤。
[0268]
在本实施例中，至少一个处理器可以构成具有对一个或多个输入执行逻辑运算的电路的任何物理设备。例如，至少一个处理器可以包括一个或多个集成电路(ic)，包括专用集成电路(asic)、微芯片、微控制器、微处理器、中央处理单元(cpu)的全部或部分、图形处理单元(gpu)、数字信号处理器(dsp)、现场可编程门阵列(fpga)或者适于执行指令或执行逻辑运算的其它电路。由至少一个处理器执行的指令可以例如被预加载到与控制器集成的或嵌入在控制器中的存储器中，或者可以存储在分离的存储器中。存储器可以包括随机存取存储器(ram)、只读存储器(rom)、硬盘、光盘、磁介质、闪存，其它永久、固定或易失性存储器，或者能够存储指令的任何其它机制。可选的是，至少一个处理器可以包括多于一个处理器。每个处理器可以具有相似的结构，或者处理器可以具有彼此电连接或断开的不同构造。例如，处理器可以是分离的电路或集成在单个电路中。当使用多于一个处理器时，处理器可以被配置为独立地或协作地操作。处理器可以以电、磁、光学、声学、机械或通过允许它们交互的其它手段来耦合。
[0269]
在本实施例中，非临时性计算机可读存储介质可以是rom、随机存取存储器(ram)、cd
‑
rom、磁带、软盘和光数据存储设备等。
[0270]
此外，本发明的特征和益处通过参考示例性实施例进行说明。相应地，本发明明确地不应局限于这些说明一些可能的非限制性特征的组合的示例性的实施例，这些特征可单独或者以特征的其它组合的形式存在。
[0271]
本领域技术人员在考虑说明书及实践这里公开的公开后，将容易想到本公开的其它实施方案。本技术旨在涵盖本发明的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本发明的一般性原理并包括本发明未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本发明的真正范围和精神以权利要求书为准。