地质系统的空间结构与动力学过程的计算模拟及其在成矿预测中的应用.doc

第5讲 地质系统的空间结构与动力学过程的计算模拟及其在成矿预测中的应用主讲:刘亮明5.1 计算模拟的重要性Because many geological phenomena and processes are based on complex regularities and natural laws which are additionally overlain by numerous random influences, temporal alterations and/or spatial deviations, to understand and illustrate them without some modeling is nearly impossible H. Thiergrtner(国际数学地质学会主席)Geological evolution and geodynamic processes cannot be adequately reproduced and studied in a laboratory, not with real rocks subjected to elevated temperatures and pressures，nor with analogue materials at more convenient conditions. Simple aspects of real rock behaviour can be produced, but not complex behaviours involving multiple deformations, metamorphism and phase changes at high temperatures and pressures, fluid flows and coupled chemical reactions. The only way that complex geological events and geodynamic processes can be studied is via numerical simulation on a computer. Such simulations require robust numerical codes that can reproduce the full range of real rock behaviours, e.g. elastic, plastic, brittle, ductile, creep, bulk and localised material flow, granular flow, folding, fracture, compaction, dilation, solid-state transformation, mineral growth and dissolution (metamorphism), melting, fluid flow, etc. Also, any strong feedback must be incorporated so that each geodynamic process (deformation, thermal flux, fluid flow and chemical reaction) can influence the operation of the other processes.G. P. Price（澳大利亚地球动力学合作研究中心主席）Computational modeling is a promising method for investigating the spatial, topographic characteristics and geodynamic processes of complex geological systems involving irregularly-shaped geological bodies.刘亮明125.1.1 计算模拟是地球科学发展的必然途径地质调查和研究的根本目的是什么？答案：4 WWhat are the geological phenomena and/or the geological bodies, including the composition and architecture?Where were they formed?When were they formed?Why do they occur in such a pattern?4W实际上就是地质系统的空间结构和动力学过程。地质研究的根本任务就是了解地质系统的空间结构和动力学过程。 地球科学的最致命的缺陷：预测能力不足。计算模拟是提高地球科学预测能力的有效途径之一，如同大气动力学计算模拟对提高天气预报的准确性一样。5.1.2地质系统的特殊性地质系统的特殊性：（1）复杂，地质要素的形态和拓扑关系的复杂，（2）不透明，不仅肉眼不可能见到每一点的情况，也没有哪种仪可以短时间内探测清每一点的情况；（3）巨大的空间尺度，但探测的数据非常有限，而且可能测点也不规则不均匀分布；（4）巨大的时间尺度，时间上多过程叠加；（5）复杂的动力因素，不同过程耦合，最常见的MTHC耦合；（6）极端的物理化学条件，高温、高压等极端条件；地质系统的这些特殊性决定：我们不可能直接观察记录所有地质系统，也不能在实验室再造一个完整的地质系统的演化过程，要了解地质系统的空间结构和动力学过程，没有模拟是不可能的。计算模拟是最佳的模拟手段这一。5.1.3 计算模拟的优点科学研究的三大方法：（1）理论分析：通过理论推理，能表示成数学形式（2）实验：获取实验数据，总结规律，传统的地球科学主要靠经验积累，也是一种实验。（3）计算模拟：用计算的方法，在计算机的虚拟时空中再现地质系统及其演化过程。通过计算模拟获取科学知识，包括规律和形成规律的原因。较之传统的地质研究方法，计算模拟具有如下的优点：（1）成平低，速度快，易重复；（2）时空过程可监控；（3）条件不受限制；（4）其理论同内核强于传统的地质科学；（5）定量，避免了传统地质定性推理的随意性。5.2 地质系统空间结构三维模拟的方法5.2.1 地质系统空间数据的特点及模拟方法的分类在三维空间内，地质系统的地质要素其实主要只有两种：面要素，包括地层界面、不整合面、断层破裂面、岩体表面等；体要素，包括岩体、地层、矿体等。我们对地质系统的了解是通过地质观察、物化探和钻探等获得的，其实主要是一些点数据，这些点数据具有如下特点：（1）数据点非常有限，大多数情况下是不充分；（2）点的分布并不是均匀规则的； 我们要了解和表征地质系统，其实就是通过这些不均匀分布、也不充分的点源数据构成面，然后再由面构成体。传统的地质研究其实只能由点构成线，然后再构成各种二维图：平面图、断面图和剖面图。然后再由这些二维的图去想象三维空间的形态。很显然这种“二维平面 + 空间想象”的表达方法具有极大的局限性，不足以高效率地准确表达复杂的地质系统。空间结构的计算模拟则是通过计算的方法将点源数据构成面，然后由面构成体。根据由点源数据生成面的方法的不同，计算建模的方法分为显式模拟(Explicit Modelling)和隐式模拟(Implicit Modelling)。5.2.2 三维构造的显式模拟（1）显式模拟的基本操作显式模拟即手工圈定边界，是目前国内常用的模拟方法。一般三维勘探软件的矿体建模和地质建模都是通过这种方法进行，由剖面图上线连成面（一系列的三角面），实际上是剖面上的点相互连接构成一系列相互连接的三角面，封闭的曲面即构成一个体，所以其基本步骤包括：（1）由点和线直接生成面，其实是由3个相临的点构成一个三角面，由一系列的三角面模拟一个曲面；（2）由封闭的曲面构成体；（3）圈闭的面构成体和体的剖分。（2）面的光滑和体的剖分在相临的数据点之间插值，将有棱角的三角面转化为光滑的曲面。GOCAD的DSI （discrete smooth interpolation）插值算法。体的剖分就是将一个体划分成一系列四面体或六面体的体积单元。（3）应用G oCAD进行显式模拟的基本操作(1) Creating the 2D outlines of the geological body from cross-sections or contour maps (图5-1a); (2) Connecting the adjacent section outlines by using a triangulated irregular network (TIN) (图5-1b); (3) Optimizing the TIN through discrete smooth interpolation (DSI) in order to construct a set of equilateral triangles (图5-1c); (4) Deleting the triangle mesh and outputting the model surface as smooth as the real geological interface (图5-1d)图5-1 应用GoCAD建立面模型的步聚5.2.3 三维属性的隐式模拟（1）隐式建模的基本操作隐式模拟即通过已有采样点数据，确定分布函数，然后通过分布函数由已知推未知，将不完备的样品数据通过空间插值生成完备的样品数据，然后再通过三维曲面构建算法自动生成三维可视模型。由于显式模拟工作量大，非常费时，且主观随意性太强，不同的人圈出来的模型经常差异很大。随着Kriging、RBF等隐式算法的发展和成熟，隐式模拟在自动生成、快速适应新数据以及结果的精确性方面则越来越显示出优越性（2）数据的空间插值三维属性隐式模拟的核心是数据的空间插值，插值的方法很多，常用的有：距离反比法（Inverse Distance Weighted，IDW）、趋势面（Trend Surface ）法、样条曲线法（Spline）、克里金法（Kriging）、径向基函数法（Radial Basis Function，RBF）。其中Kriging法是地质建模中最常用的插值方法。Kriging法是由南非采矿工程师Kriging在研究金矿时首次提出，然后由法国数学家马特龙加以发展和完善。它是一种无偏的、误差最小的、最优化的主要用于地学领域的基于统计理论的方法。其基本原理是，任一待估点V的真值的估计值是估计邻域内n个信息值的线性组合，即 (1)式中待估邻域内每个样品权值i的求解方程组为:K=M (2)式(2)中的K，M矩阵如式(3)： (3)式(3)中的为拉格朗日系数，为i，j两点的变异函数值，求解方程为：h=0(h)= 0h式4)中的0为块金常数(Nugget)，为基台值(Sill) ，为变程(Range)， h是i，j两点间的空间距离。最简单的情况如：变异椭球体和X，Y，Z轴重合且各向同性，即变异椭球体三轴变程相同(为球体)，则h = (x2-x1)2+(y2-y1)2+(z2-z1)2 (5)如变异椭球体和X、Y、Z轴不重合，则在求距离前，需先对(x,y,z)坐标做旋转变换，使之和变异椭球体重合,当变异椭球体各向异性，即X、Y、Z三个方向上的变程不同时，还需对变异椭球体进行变换，使之变为各向同性，然后再计算两点间的空间距离。求出两点之间的空间距离h后，代入式(3)，即可得到两点之间的变异函数值，再将之代入式(2-2)，即可求出邻域内每个样品的权值i，再将求得的i代入式(1)，即可求得待估点品位值。对每个待估点依次重复上述步骤，即可求得每个待估点的估计品位值。（3）等值面的提取从插值后的数据中提取等值面是隐式模拟的关键步骤之，生成等值面的算法很多，其中Marching cube算法是最常用的一种。Marching cube算法是三维数据场等值面生成的经典算法，Marching cube算法的基本思想是逐个处理数据场中的立方体，分类出与等值面相交的立方体，然后再通过插值计算出等值面与立方体的交点。Marching cube算法首先对立方体的八个顶点进行分类，以判定顶点是位于等值面之外，还是位于等值面之内，再根据八个顶点的状态，确定等值面的模式。然后直接由立方体各顶点的状态检索出等值面的分布，确定该立方体体素内的等值面三角片连接方式。在确定立方体内三角片等值面分布模式后，即可用线性插值计算等值面与立方体边的交点X:X=i+(Value-Fi,j,k)/(Fi+1,j,k-Fi,j,k) (6)式(6)中的Value为阈值，Fi,j,k为顶点(i,j,k)的数据值由于函数F(x,y,z)的梯度垂直于等值面，因而可用数据场Fi,j,k在等值点P的梯度做为P点处的法向量，又因数据场是三维正交数据场，故可以采用中心差分计算出立方体顶点的梯度，再用线性插值计算P点的梯度值: (7)式(7)中的，为立方体的大小，为计算得到的顶点(,)的梯度值。5.2.4 空间结构模拟的实例（1）浅成岩体及相关矿体的三维形态模拟根据钻探和地球物理测深资料模拟某岩体和相关矿体的三维形态，揭示它们之间的空间关系（图5-2）。图5-2 某岩体及相关矿体的三维模拟（2）斑岩型铜钼矿的品位的三维空间分布模拟根据钻孔分析数据，模拟某斑岩矿化体系中Cu和Mo三维空间分布，揭示Cu和Mo的空间聚集规律，进行矿体边界品位优化（图5-3）。5.3 地质系统动力学过程的计算模拟5.3.1 动力学计算模拟的理论基础动力学的核心数学问题是动力学方程组，主要由守恒方程(能量、质量和动量)和系统的本构方程组成，实际上是一组复杂的偏微分方程组成，理论上说研究这些方程组应该有两种不同的方法，一是解析法，或称分析法（Analytic Method），另一种是数值法（Numerical Method）或称计算方法（Computational Method）; 解析方法就是用数学上的积分法或积分变换等方法直接求数学模型的解，其解称为解析解，它是数学模型的精确解，这种解的最大优点在于把表征物质运动规律的变量与激发条件、时空变化包含在一个表达式中，这样便于分析物质运动规律。但实际上由于动力学方程的复杂性和系统几何结构和形态的复杂性，大多数情况下，我们是不可能求得分析解的；Mo 0.02Mo 0.03Mo 0.04Mo 0.05图5-3 某斑岩铜钼矿床Mo的空间分布的三维模拟计算模拟数（Computational Modeling）又称数值模拟(Numerical Modeling)就是离散方法求解数学模型，或者说用数值方法求解动力学方程组的数值解，其解为数值集合，是数学模型的近似解。尽管不如解析法精确，但它对复杂的动力学系统而言，是一种切实可行的方法。实际上对地质研究而言，其精度已足够高；计算模拟实际上就是在计算机的虚拟空间内定量地再造某一自然事件或回放真实的过程(quantitative reconstruction of a natural event or a replay of a “real-life” process)。描述地质系统动力学的常用的数学方程： (8) (9) (10) (11) (12) (13) (14)5.3.2 动力学计算模拟和基本步骤进行地质系统的动力学计算模拟一般都按如下7个步骤进行：1）提出问题 A “Story” or key question；2）软件 A numerical code to simulate this story3) 构建几何模型 Construction of geometry (mesh)4) 本构关系 Selection of appropriate rheology (constitutive laws)5) 模型组成单元的特征参数 Specification of relevant rock properties (parameters) for your model；6) 边界条件和初始条件 Defining necessary boundary conditions for the model reflection of natural deformation or thermal environment；7) 模型运算及结果输出，Run the model computational iteration。5.3.3 动力学计算模拟的难点1）算法与合适的软件，没有完全适合于地质系统的多过程耦合动力学计算模拟软件；2）复杂系统的有限单元模型的建立，建立复杂地质系统的六面体单元模型非常难；3）单元参数与初始及边界条件的确定。5.3.4 动力学计算模拟的算法与软件（1）动力学计算模拟的算法常用的算法有：（1）有限元(Limited Element)法；（2）有限差分(Limited Difference)法；（3）边界元 (Boundary Element) 法；（4）离散元（Discrete Element）法等。不管是哪种方法，都得将一个连续的系统分割成n单元，如果整个系统有m个独立的变量的话，那么每算一步，就要解mn个方程。如果系统尺度很大，而我们要保证模拟的准确度，那么n就必须足够大，这样就导致计算工作量相当大，在计算机技术不发达时，这对数值模拟来说是一个相当大的难题。就现代计算机技术水平及发展趋势来看，这已经不是一个大的问题。（2）动力学计算模拟的软件复杂地质系统的动力学计算模拟一般难以依靠研究者独自开发软件来完成，一般是依靠一些大型商业软件来进行，在这些商一软件平台的基础上进行二次开发。常用动力学计算模拟软件包括：（1）FLAC、PFC和UDEC (/), （2）ANSYS (/)；（3）FLUENT (/)；（4）COMOSOL (/)等。FLAC（Fast Lagrangian Analysis of Continua，连续介质的快速拉格朗日法分析）是我们常用到的一种软件，由美国Itasca公司开发的一种以显式有限差分法（Explicit Finite Difference）进行连续介质物理作用模拟的商业软件，它源自数学力学的拉格朗日元法。最初由Minnesota大学的Peter Cundall开发的。目前，FLAC已由二维发展到三维，二维计算程序V3.0以前的为DOS版本，V2.5版本仅仅能够使用计算机的基本内存（64K），所以，程序求解的最大结点数仅限于2000个以内。1995年，FLAC2D已升级为V3.3的版本，其程序能够使用护展内存，大大发护展了计算规模。FLAC的基本原理和算法与离散元相似，但它却像有限元那样适用于多种材料的模拟和边界条件非规则区域的连续问题求解。在求解过程中，FLAC采用了离散元的动态松弛法，不需要求解大型联立方程组（无需形成刚度矩阵）；另一方面，同以往的差分分析法相比，FLAC不但可以对连续介质进行大变形分析，而且还可以模拟岩土体沿某一软弱面产生的滑动变形；FLAC还能针对不同的材料，实现相应的本构方程来比较真实地反映实际材料的动态行为。程序将计算区域内的介质划分为若干个二维单元，单元之间用节点相互连接。上述过程中，网格的划分与有限元不同之处在于其网格分物理网格和数学网格且互为影射，所划分的网格只要有序也可具有不规则的形状。5.3.5动力学计算模拟的实例以FLAC3D为平台，模拟冬瓜山矿田的青山脚岩体在同构造冷却过程的力热流耦合动力学演化，揭示了冬瓜山铜矿的层状形成是在岩体的同构造冷却过程中变形和流体作用的结果（图5-4）。图5-4 青山脚岩体同构造冷却MTH耦合结果5.4 计算模拟在成矿预测中的应