模型分析

1. 模型分析

由于影响地热流体动态变化的因素非常复杂,具体模型已经过高度概化,与实际情况有一定误差,然而使模型参数无限接近含水层固有参数,也是不必要和不可能的。产生差异的影响因素可能有:①某些观测孔地热流体水位不能准确代表该层实际值。特别是地热流体温度和密度对水头高的影响,虽然对地热流体作了统一温度等技术处理,但是否合理还有待在实践中进一步验证;②尽管十分注重长观孔的观测质量,难免存在观测误差;③长观及统测数据存在偶然性,不能代表某一时段的水位;④因剖分的关系,计算水位标高代表的是整个单元的平均水头高,而观测值是某个点的值;⑤实测地热流体水位等值线图是人工绘制的流场,存在人为误差;⑥地热地质模型概化产生的误差及模拟过程的误差。
数值模型存在来自多方面的不确定因素。为了使所建模型确实能代表所研究的地质体、再现野外实测水位和流场,必须对所建模型进行分析评估,根据建模的目的要求,对模型变量之间的依赖关系、稳定性、系统参数的灵敏度等进行分析。通过分析、确定不确定因素对校正模型的影响程度。如果不符合要求,就修改或增减建模假设条件,重新建模,直到符合要求。
(1)误差分析
实测热流体水头和模拟水头等值线图的对比提供了一种看得见的定量方法,同时给出了识别误差空间分布的大致信息。但根据野外资料绘制的等值线图包含了绘制引起的误差,因此不宜作为识别的唯一证据。误差是衡量实测值和模拟值之间接近程度的物理量。实测值与计算值之间的误差大小和某种误差值的平均是展示识别结果的常用方法。
A.平均误差(绝对误差)
误差值的平均作为识别过程中的平均误差,识别的目的就是使这个误差成为最小。本次工作采用均方根(RME)误差和平均绝对误差(MAE)来衡量水头实测值和模拟值间的接近程度。其计算公式如下:

沉积盆地型地热田勘查开发与利用

式中:i为表示比较时刻(第i次实测);n为观测井实测次数;hob为实测热流体水头(m);h为计算热流体水头(m);RME为均方根误差值;MAE为平均绝对误差值。
B.平均相对误差
相对误差体现了误差值对整个系统内总水头变化的影响程度,一般用它的百分比表示,其计算公式为

沉积盆地型地热田勘查开发与利用

式中:Er为平均相对误差;i为表示比较时刻(第i次实测);n为观测井实测次数;hob为实测热流体水头(m);h为计算热流体水头(m)。
由上式可看出Er值越小,说明模型拟合得越好;反之Er值越大,说明误差值对整个系统内总水头变化的影响程度越大。
本次工作利用公式5-26至公式5-28对模型模拟期和检验期部分长观井实测值和计算值间误差进行分析计算。但数值模拟无法给出一个绝对指标来衡量计算结果的好坏。从误差分析结果来看,模型存在一定的误差,但都不算大。根据其总体误差分析值判断,模拟值与实测值拟合较好,模型模拟基本达到精度要求。
(2)灵敏度分析
由于水文地质参数、储层参数、岩性、构造、外部影响以及边界条件都不可能知道的很详细,因此建模过程中通过判断所取得的上述参数存在着不确定性。灵敏度分析的目的是为了对已经识别过的模型的不确定性进行量化,用来考察微观变化对建立模型的整体影响,检验模型的优劣性。
灵敏度分析包括局部灵敏度分析和全局灵敏度分析(包括定量全局灵敏度分析、定性全局灵敏度分析),本次工作采用国内使用较多的局部灵敏度分析方法。局部灵敏度分析也称一次变化法,其特点是在其他参数不变的情况下,只针对某一个参数,评价模型结果在该参数每次发生变化时的变化量。通常采用灵敏度系数作为衡量参数灵敏度的标准,其形式是:

沉积盆地型地热田勘查开发与利用

式中:βi,k为模型的因变量(如水头、浓度等)对第i观测点,第k参数的灵敏度系数;Hi为模型第i观测点的因变量(此处指热流体水头);ak为模型第k参数(可以是模型中用到的任一参数)。
由于不同模型变量的量纲不同,不同参数的量纲差别也很大,所以式5-29所示灵敏度没有一个统一的单位和量纲。为了便于不同参数间灵敏度的比较,式5-29也可化为下列形式:

沉积盆地型地热田勘查开发与利用

这样灵敏度的量纲就和模型变量的量纲一致(式5-30中参数参见式5-29)。
βi,k值大意味着第k个参数值改变对模型变量(热流体水头)的影响大,针对某一特定参数k的灵敏度系数值,在保证其他所有参数不变时,使该参数k的值由ak变化为ak+Δak相应的因变量值由Hi(ak)变化为Hi(ak+Δak),由式(7 11)获得其近似值,即:

沉积盆地型地热田勘查开发与利用

或采用下列标准化的形式:

沉积盆地型地热田勘查开发与利用

式中:Δak为模型第k个参数的细小改变,Hi(ak)和Hi(ak+Δak)表示第k个参数值为ak和ak+Δak时第i观测点所得的模型变量值(热流体水头)。
式5-30和式5-32给出的灵敏度只是在一个特定位置上对一个给定的参数,模型所反应的灵敏度的大小。
目标函数均方根(RMS)误差对某一个模型输入参数的灵敏度定义了一个单一的灵敏度,因为它不会由于模型变量不同而有不同的灵敏度。为此以目标函数的均方根(RMS)误差来代替灵敏度公式中的模型变量,有:

沉积盆地型地热田勘查开发与利用

或标准化形式:

沉积盆地型地热田勘查开发与利用

式5-33意味着βk≈  ,其中ΔH为由于改变参数值ak为ak+Δak时,目标函数均方根(RMS)误差。
选取模型中不确定因素较多的水文地质参数:水平向渗透系数Kxy、垂向渗透系数Kz和弹性释水系数S*分别增加和减小10%,20%,利用公式5-34来计算目标函数的灵敏度和对数学模型进行灵敏度分析。由于模型中参数分区赋值,灵敏度计算时ak和ak+Δak取其平均值,计算结果见表5-8。

表5-8 Ng热储层Kxy,Kz,S*灵敏度分析计算表

由上表可知,馆陶组灵敏度系数最大的是弹性释水系数的2.50;水平向渗透系数和垂向渗透系数的灵敏度系数分别为1.72和1.71,基本相当。
从灵敏度系数值的大小可以看出,各参数的取值变化对计算结果有一定的影响。但影响程度较小。可见在正确选择各项水文地质参数之后,均表现为不太灵敏的参数。
进行灵敏度分析是为了确定模型的结果对模型参数的敏感程度,如果模拟结果对某一特定参数高度敏感,模型作出重要解释和预报的能力将受到和该参数有关的不确定性的严重影响。根据表5 8的结果,说明所建模型相对比较稳定,模拟结果对给定参数不太敏感,也即意味着水文地质参数领域的不确定性对模型解释和预报能力的影响有限。
(3)流量均衡分析
从均衡的角度出发,热流体均衡是对计算结果的可信度的一个重要测量指标。有限差分方程是依据流体的连续性方程建立起来的,所以流进和流出一个地下水系统的流量总和应能满足连续性原则:总流入量和总流出量之差等于贮水量的变化量。模拟的热流体均衡变化与实际要基本相符。
A.均衡区的确定及均衡时段划分
均衡区的范围以明化镇组和馆陶组计算范围为界,均衡计算采用多年平均及2002～2006年来进行均衡计算评价。
B.均衡方程的建立
地热流体资源一般埋藏较深,相应补径排项较为单一。依据均衡原理,结合热流体补给、径流、排泄条件,建立均衡方程如下:

沉积盆地型地热田勘查开发与利用

式中:Q补为地下水总补给量(m3);Q排为地下水总排泄量(m3);S*为弹性释水系数;F为均衡区面积(km2);Δt为均衡时间段长(a);ΔH为与Δt对应的水位变幅(m);Q径入为侧向径流补给量(m3);Q开采为人工开采总量(m3);Q径出为侧向径流排泄量(m3);Q越入为此处主要指热储层通过导水断裂得到的越流补给量(m3);Q越出为热储层越流流出量(m3);Q泥岩释水为热储层中泥岩层压缩变形释放出的水量(m3);Q回灌为人工回灌总量(m3)。
C.热流体均衡分析
通过上述热流体各均衡要素分析可知评价区补给、排泄项均比较单一,地热流体系统模拟期热流体量均衡结果见表5-9。

表5-9 Ng组热储层模拟期均衡表

从表5-9中可以得到,Ng模拟期总补给量为2913.47×104m3,包括侧向补给量、弹性释水量、泥岩释水量和部分从沧东断裂带得到的垂向越流补给,分别占总补给量的54.85%,23.52%,4.87%和16.81%。弹性释水量动用的是储存量,因此弹性释水量越多,水位下降越大。侧向补给量在模型中是通用水头边界控制的流入流出量。从各补给项比例关系可以看出,评价区主要的补给来源是侧向补给。
通过上述分析可以看出,水均衡计算仅反映模型计算过程的可靠性,并不能反映模型与野外实际条件的一致性。

2. 模型研究分析

薄层的差异地震响应呈现为简单波形，厚的薄互层的差异地震响应则表现为复合波形，二者对应的 差异波形有明显的不同，故需要分开来分析。

图7.6 薄层基本差异地震模型的速度变化解释 a—增速模型；b—减速模型；c—速度不变模型

（1）薄层差异波形解释
薄层的差异地震响应是简单的差异波形，可以对差异波形进行简单解释。通过正演模拟得知，简单 差异波形与薄层中砂岩和围岩的波速大小无关（见图7.1），只与砂岩的波速变化有关，因此，将薄的薄 互层的对应的差异波形大体分为三类来分析，即正波形、负波形和零波形，建立三种差异解释模型（图7.6），差异数据体中的正波形对应砂岩波速增大，负波形对应砂岩波速减小，零波形对应砂岩波速没有变化。总厚度下与四分之一波长的薄互层的差异波形应与薄层的差异波形相似。

图7.7 厚的薄互层对应的复杂波形

（2）厚的薄互层差异波形解释
由于厚的薄互层的地震差异响应是多峰复合波（图7.7），波形的组合比较复杂，对应的波速变化组 合也不是确定的，因此，在对差异剖面中的多峰复合波进行解释时应该谨慎，只有依据地质、测井资料 建立模型，对各层的波速变化进行正演模拟，以此来解释复杂的多峰复合波。
（3）厚的薄互层条件下顶部油层的差异波形具有显著屏蔽效应
通过模型研究认识到，顶部油层的差异响应可以掩盖下部油层的差异响应。如图7.8，以五层的砂泥 岩互层来模拟厚的砂泥岩薄互层，泥岩是高速地层，而砂层是低速油层，图7.8（a）是顶部砂岩增速、底 部砂岩减速模型，图7.8（b）是（a）模型的差异地震响应，分析楔形模型的波形响应，可以分为三种情况：

图7.8 屏蔽效应模型（一）

第一种情况，当中间的泥岩夹层厚度远远大于调谐厚度时，顶部砂岩增速对应完整的正波形，而底 部砂岩减速，对应完整的负波形，并且顶部与底部砂岩的差异地震响应完全分开。这种情况下，通过差 异波形可以判断顶、底部油层的波速变化。
第二种情况，当中间的夹层厚度小于调谐厚度时，首先，顶部砂岩增速对应的正波形与底部砂岩减 速对应的负波形的旁瓣波谷互相叠合形成多峰复合波，这种情况下，通过差异波形正演模拟可以判断顶、底部油层的波速变化。
第三种情况，顶部砂岩增速对应的正波形与底部砂岩减速对应的负波形完全叠合形成单峰复合波。这种情况下，顶部砂岩的波速变化完全屏蔽了底部砂岩的波速变化，通过差异波形只能判断顶部砂岩的 波速变化，而不能识别出底部砂岩波速的减小。
如图7.9，以五层的砂泥岩互层来模拟厚的砂泥岩薄互层，泥岩是高速地层，而砂层是低速油层，图7.9（a）是顶部砂岩减速、底部砂岩增速模型，图7.9（b）是（a）模型的差异地震响应，分析楔形模型的波 形响应，可以分为三种情况：

图7.9 屏蔽效应模型（二）

第一种情况，当中间的泥岩夹层厚度远远大于调谐厚度时，顶部砂岩减速对应完整的负波形，而底 部砂岩增速，对应完整的正波形，并且顶部与底部砂岩的差异地震响应完全分开。这种情况下，通过差 异波形可以判断顶、底部油层的波速变化。
第二种情况，当中间的夹层厚度小于调谐厚度时，首先，顶部砂岩减速对应的负波形与底部砂岩增 速对应的正波形的旁瓣波谷互相叠合形成多峰复合波，这种情况下，通过差异波形正演模拟可以判断顶、底部油层的波速变化。
第三种情况，顶部砂岩减速对应的负波形与底部砂岩增速对应的正波形完全叠合形成单峰复合波。这种情况下，顶部砂岩的波速变化完全屏蔽了底部砂岩的波速变化，通过差异波形只能判断顶部砂岩的 速度变化，而不能识别出底部砂岩波速的减小。

3. 计算模型的建立

5.3.3.1 地质条件的概化
为方便起见，坐标设定基岩面为零点。对塌陷体进行地质概化后，基本可分为如下三层：
（1）土层及耕植土层，2.1～4.0m。
（2）砂质粘土层，0～2.1m。
（2）基岩层，z=0～18m，应属于阻-透型地质概化模型。
5.3.3.2 几何模型及边界条件
几何模型分为两种情况进行模拟。第一种情况，模型没有抽水井，而是在相应的边界上加上相应的流量及降深，大小为40m×30m×18m，以边界排水来进行抽水井的模拟，模型外观见图5-9。另一种情况，则模拟了抽水井的存在，目的在于观测整个水面降落漏斗的变化。大小变为60m×30m×18m，按实际抽水量换算后，模拟流量为0.023m3/s模型外观见图5-10。边界条件设计为进水边界，两种情况下的边界流量在下面的补给边界中列出。设计岩溶洞穴内有水压及水流存在。坐标原点定在基岩面上，向下为负。边界情况如下。
静力约束边界：
约束x=9.9～10.1m所决定的面，使之在x方向位移为零。
约束z=-18.1～-17.9m所决定的面，使之在z方向位移为零。
约束x=69.9～70.1m所决定的面，使之在x方向位移为零。
约束y=-14.9～-15.1m所决定的面，使之在y方向位移为零。
约束y=14.9-15.1m所决定的面，使之在y方向位移为零。
补给边界：
井壁：0.03m3/s
流量边界-1×10-5m3/s·m2x=69.9，70.1 y=-15，15 z=-18，0
流量边界-1×10-5m3/s·m2x=9.9，10.1 y=-15，15 z=-18，-15
流量边界-1×10-5m3/s·m2x=10，70 y=14.9，15.1 z=-18，3
流量边界-1×10-5m3/s·m2x=10，70 y=-14.9，-15.1 z=-18，3
流量边界-1×10-5m3/s·m2x=10，70 y=-15，15 z=-18.1，-17.9
流量边界0m3/s·m2 x=69.9，70.1 y=-15，15 z=3，4.1
流量边界0m3/s·m2 x=-9.9，-10.1 y=-15，15 z=3，4.1
流量边界0m3/s·m2 x=10，70 y=14.9，15.1 z=3，4.1
流量边界0m3/s·m2 x=10，70 y=-14.9，-15.1 z=3，4.1
数值模拟中，在土层与基岩接触的地方设计了一个土洞，直径1m。在土洞的下部设计了一个岩溶洞穴，以模拟基岩中岩溶的发育。试验结果表明，该岩溶洞穴在水位下降致塌中起着重要作用。模型外观图见图5-9及图5-10。实现数值模拟的数值语句见附录3。

图5-9 无井边界排水模拟外观图


图5-10 有井排水模拟的模型外观图

5.3.3.3 参数的选取及计算模型
FLAC3D可以解决一般的水流力学问题。在模拟水流力学问题中，FLAC3D具有以下特征：
（1）可计算在地下均匀渗流的力学过程，流域内的不透水材料被命名为空模型（null model）；
（2）不同的单元群可以被命名为不同的流动模型或属性；
（3）可以通过命令定义补给及不透水边界，或漏水边界；
（4）地下的集中排泄或是补给可以通过定义模型中的节点或区域来实现；
（5）计算可以有显式计算或隐式计算两种。
FLAC3D中的数学模型以Biot理论为基础，可以描述不同类型的流动过程，包括空隙介质中的达西流（Darcy），以及空气及水的流动。用来定义流量及孔压之间关系的表达式为达西定律：
qi=-k［p-ρfxjgj］i
式中：k为渗透系数（或流动系数）（m4/N·s）；ρf为流体密度（kg/m3）；gi，i=1～3，为重力矢量的三个元件（m/s2）。
FLAC3D中使用的特殊参数Biot系数及Biot模量与土元件的体积模量有关，其表达式为：

岩溶塌陷机理及其预测与评价研究

式中：K为材料的排水体积模量，Ks为土元件的体积模量，其大小为：  到1之间。
Biot模量为：

岩溶塌陷机理及其预测与评价研究

式中：Ku为材料的不排水体积模量。在FLAC3D中还需要输水的体积模量，对于非压缩材料，有：

岩溶塌陷机理及其预测与评价研究

式中：Δp为水压的变化值；ΔV/V为体应变。而对于水的压缩性，在室温下Kf=2×109Pa。FLAC3D可以对细粘土中所产生的拉力加以限制。在很细的粘土中，水加入其中时，可以产生拉力，在FLAC3D中用负压表示，但这并不同于毛细水产生的负孔压。这是由于水的充填使材料的体积膨胀而引起，有相应的命令可以去掉这个负压。
对于实际问题，参数的选取主要根据实际的岩性及试验数据，参考《工程地质手册》而选取。参数包括两部分，第一部分为静力计算所需的材料性质参数；第二部分为涉及水流的力学参数，包括渗透系数及孔隙度等。对于后者其大小将直接影响水流在介质中的运动性质。模型中各层的第一部分参数取值见表5-1。

表5-1 静力材料参数

所用的力学计算模型，对于建筑层及基岩采用弹性模型（elastic），对于硬塑、可塑粘土层、砂质粘土使用摩尔-库仑模型（Mohr-Coulomb）（钱家欢等，1991；雷晓燕，1999）。
模型取值如下。
重力加速度：-9.8m/s2。
重度：
水，1000kg/m3
粘土，2400kg/m3，z=2.1～2.9m。
砂质粘土，2300kg/m3，z=0～2.1m。
基岩：2700kg/m3，z=-18～0m。
粘土：2100kg/m3，z=2.9～4m。

4. 计算模型建立

7.4.1.1 29210工作面(剖面1)
模型取宽度为700m，高度为340m。模型共有33000个平面单元，网格尺寸平均为3m×3m。模型采用位移边界条件，即模型的左侧、右侧边界水平位移固定，地面垂直位移固定。

图7-59 29210剖面模型

7.4.1.2 29211工作面(剖面2)
模型取宽度为1000m，高度为240m。模型共有24000个平面单元，网格尺寸平均为3m×3m。工作面计算模拟模型如图7-60。模型采用位移边界条件，即模型的左侧、右侧边界水平位移固定，地面垂直位移固定。
7.4.1.3 北帮三维模型
模型走向长度为850m，倾向长度为700m，高度为340m，共剖分166999个立体单元，模型开采4#和9#两层煤，工作面推进方向沿着煤层的倾向方向，即“顺向开采”。模型采用位移边界条件，即在三维数值分析中，模型的东西南北四个铅垂面的水平位移固定，地面垂直位移固定。模型的地质分层情况与剖面1、2相同。岩土体的本构模型采用Mohr-Coulomb模型，对采空区采用开挖模型(null)。

图7-60 29211剖面模型


图7-61 安家岭露天矿北帮三维分析模型

5. 模型建立及预测

根据对影响各煤层甲烷含量地质因素的分析研究，我们从中选出了三至四个因素指标作为影响甲烷含量大小的主要控制因素，作为建模预测的指标（表8.2），根据GM（0，N）建模原理，利用已知瓦斯钻孔获得的各煤层吨煤甲烷含量值与相应的主控因素指标统计值（表8.3至表8.8），便可建立起GM（0，N）预测模型：

表8.2 韩城矿区各煤层甲烷含量主要影响因素一览表

北区：
2＃煤层为：
Qd＝1.238998＋1.849812×10－2MCMS－0.4038353DBYX－0.3008947T20M
3＃煤层为：
Qd＝5.284516＋0.9561405GZ－3.103036×10－2SYZH＋
3.650546×10－3MCMS－1.065613XPXS
11＃煤层为：
Qd＝12.35199－1.592682GZ－0.5038853T10M＋3.1008224×10－2DCZH
南区：
3＃煤层为：
Qd＝8.57431＋1.159906×10－3MCMS－0.1627932D10M－0.5087085T10M
5＃煤层为：
Qd＝16.54584－1.378349GZ－0.151046SYZH＋4.48066×10－3MCMS－05600767D5M
11＃煤层为：
Qd＝8.730469＋7.424534×10－3MCMS－1.633005DBYX－0.5197687D10M

表8.3 北区2号煤层甲烷含量与主控因素指标统计表


表8.4 北区3号煤层甲烷含量与主控因素指标统计表


续表


表8.5 北区11号煤层甲烷含量与主控因素指标统计表


续表


表8.6 南区3号煤层甲烷含量与主控因素指标统计表


表8.7 南区5号煤层甲烷含量与主控因素指标统计表


续表


表8.8 南区11号煤层甲烷含量与主控因素指标统计表

利用所建立的预测模型，分别对各煤层已知钻孔甲烷含量值进行了预测，并以此计算了各孔煤层气含量实际值与预测值之间残差及平均相对误差值（表8.9至表8.14）。从表中看出，尽管各煤层钻孔相对误差值较高，但总体相对误差平均值均＜25%，其中北区2＃煤层相对误差平均值为13.65%，精度良好（二级）；3＃煤层相对误差平均值18.35%，精度较好（三级）；11＃煤层相对误差平均值19.18%，精度较好（三级）。南区3＃煤层相对误差平均值23.03%，精度合格（四级）；5＃煤层相对误差平均值22.01%，精度合格（四级）；11＃煤层相对误差平均值15.0%，精度较好（三级）。
因此，所建立的甲烷含量预测模型，可以用于南北区对各煤层甲烷含量值进行预测。

表8.9 北区2号煤层预测甲烷含量值误差检验表


表8.10 北区3号煤层预测甲烷含量值误差检验表


续表


表8.11 北区11号煤层预测甲烷含量值误差检验表


表8.12 南区3号煤层预测甲烷含量值误差检验表


表8.13 南区5号煤层预测甲烷含量值误差检验表


表8.14 南区11号煤层预测甲烷含量值误差检验表

6. 计算模型的建立

6.6.1.1 实例选取及地质条件的概化
实例的选取仍以形态规则、过程典型为原则。通过比较最后选取了位于徐州矿务局青山泉—权台煤矿专用铁路K4+315m至K4+495m地段发生的岩溶致塌作为数值模拟及数值试验的对象。
该岩溶塌陷现象发生在青权线龙须河段，该段近南北向延伸，塌陷发生在铁路的西侧路堤中，时间为1991～1992年汛期，典型的塌陷剖面如图6-2所示，基岩面位于28.5m标高处，在人工河床底部有落水洞，并与石灰岩中岩溶管道相通。基岩面以上为第四系冲积层，平均厚5m，主要由粉质粘土组成，底部透水性较好，地下水面位于基岩面附近，这为土洞的形成创造了条件。
根据研究区地质情况，可将其概化为如下几层：
（1）铁路路基，32.5～35.3m；
（2）第四系冲积层，28.7～32.5m；
（3）基岩，＜28.5m，其中溶洞发育。
6.6.1.2 参数及波形选取
6.6.1.2.1 参数的选取
数值模拟所选参数主要根据《工程地质手册》选取，所需的参数如下。
（1）材料性质参数：体积弹性模量、剪切模量、内聚力、内摩擦角、抗张强度、衰减系数（表6-1）。对于衰减系数，主要采取了瑞利衰减（Rayleigh damping）及当地衰减（Local damping）。瑞利衰减通过衰减土层的刚度去实现振波在运动过程中的衰减，它与土层系统的中心频率有关（一般为22.8 Hz），在计算中使用的瑞利衰减率为0.05～0.2。当地衰减是通过对土层单元内的质量及速度的衰减实现的，它与频率无关，衰减值一般较小，主要用于静力平衡问题，应用于动力问题较少。

表6-1 模型中各层的参数选取

（2）力学参数：参照铁道设计院所编《铁路设计手册》，取货车运行时的加速度值为0.4～0.6倍重力加速度，模拟所取值为0.6倍重力加速度。据铁道部第四设计院（1984）在《高速铁路》一书中的研究，当列车以200km/h 运行时，在路轨处产生的加速度为0.5～20m/s2，考虑到伸缩缝的影响取6m/s2较为合理，相当于重力加速度9.81m/s2的0.6倍；在路基的坡角加速度一般降至路轨处的1/10左右。
6.6.1.2.2 振动波形的选取
对于简谐振动中的有阻尼受迫振动，可用大学物理学中的公式表达（张三慧，1991）：

岩溶塌陷机理及其预测与评价研究


岩溶塌陷机理及其预测与评价研究

据土动力学原理有：

岩溶塌陷机理及其预测与评价研究

式中：β——衰减系数；ω0——系统的角频率；ω1，ω——振动频率，ε、θ——相位差；A0、A、h——振幅大小。
式（6-11）中的前一部分可称为衰减项，在振动开始后很快就衰减掉。故可以省略，维持振动的主要是后一项，上式可以简化为：
x=Acos（ω1t-θ） （6-14）
这是一明显的余弦波，从野外的直接观测中发现，铁路路轨处的最大位移一般小于1cm，因此使用0.01m作为振动在输入端的最大位移。
则：
x=0.01cos（ω1t-θ） （6-15）
θ为初始相差，它与土层系统的固有周期有关。缪忠灵等（1995）曾描述过系统固有周期，对于岩溶区，可取0.1～0.8。
经简单的换算后，取ω0=0.6。β衰减系数，取0.05～0.2之间。列车振源处的振动频率ω1=2πƒ，ƒ为1，通过实测振动次数除以振动时间而得，从而式（6-15）可以表达为：

岩溶塌陷机理及其预测与评价研究

如前章所述，实际观测到的振波波形为脉冲波形，并考虑了实际中振动波应由0开始振动，最后将上式适当变形后得：

岩溶塌陷机理及其预测与评价研究

计算中产生了较好的脉冲波，如图6-7所示。另外两种可能的波形为正弦波和余弦波，分别通过变形来实现。其数学表达式如下：
正弦波（图6-8）：

岩溶塌陷机理及其预测与评价研究

余弦波（图6-9）：

岩溶塌陷机理及其预测与评价研究


图6-7 数值模拟中的输入波形（脉冲波）


图6-8 数值模拟中的输入波形（正弦波）


图6-9 数值模拟中的输入波形（余弦波）

6.6.1.3 动力加载范围
为实现对列车振动效应的模拟，振波以6m/s2的加速度被加到模型中的铁路路基上，相当于实际中的铁路伸缩缝处，其范围如图6-10所示。

图6-10 动力模拟模型外观图

6.6.1.4 几何模型及边界条件
几何模型的建模原则主要考虑以下因素：
（1）形状尽量规则；
（2）单元数量尽量合理，以减少计算时间；
（3）边界离振源尽量远，以接近实际情况，减少反射波的影响。
在模型中设计了相应的土洞，土洞的洞壁是圆滑的，在土洞下以一圆柱状的洞穴模拟了基岩中岩溶管道的存在。模型整体如图6-10 所示，模型数值语句见附录2。坐标原点设在图中左侧边界的底部正中央处，向上为正，全区范围 24.8m×30m×11.3m。
边界条件设定如下：
为了模拟的真实性，在x=0、x=24.8m、y=-15m、y=15m、z=0 所决定的面上，设计为静止边界和吸收边界，以使振波不致因边界处产生的反射波而造成对计算结果的影响。地表为自由边界。

7. 计算模型的建立

4.4.2.1 地质条件的概化
塌陷区地质结构为典型岩溶体-土层盖层结构，属单一阻水型盖层模型。土层厚9.5m，类似的土洞塌陷在研究区有几处，将土层概化后基本可分为以下几层：
（1）人工层，有杂填土及人工填土。
（2）硬塑红粘土，在地表以下0～3.5m，局部有淤泥存在。
（3）可塑亚粘土，有植物碎片，局部淤泥。
（4）灰岩，三叠系薄至中厚层白云岩、石灰岩，溶隙、溶洞发育。
各层岩土力学性质见表4-1。
塌陷直接诱发原因为距塌陷点250m处的爆破振动，塌陷发生在6月份，在研究区可见地表水的下渗及灌入现象。渗入地下的水来源主要有：水渠、生活用水、雨水。地下水位在基岩面以下。

表4-1 各层岩土力学参数表

4.4.2.2 模拟范围及边界条件
由于主要研究对象为土层中的土洞，因此，根据圣维南原理，取边界大于3倍以上洞径。实际模拟面积为30m×30m，土洞基本位于中间。垂直方向定为坐标z轴方向，为取值方便定为向下为正，坐标零点在地面。底部边界位于地下11.5m、基岩面以下2m处。在数值模拟中，主要使用了位移约束边界，约束四周及底部的位移值（位移为零）。地面为自由边界（图4-3），边界条件如下。
位移约束边界，约束x=-28.56m所决定的面，使之位移为零；
位移约束边界，约束x=5m所决定的面，使之位移为零；
位移约束边界，约束z=11.5m所决定的面，使之位移为零；
自由边界，范围z=0m及斜坡、挡墙的地表面；
位移约束边界，约束y=-10m所决定的面，使之位移为零；
位移约束边界，约束y=10m所决定的面，使之位移为零。
4.4.2.3 几何模型
鉴于FLAC3D的功能特点，研究中采用此软件进行数值模拟研究。建模时除对研究块体进行模拟外，还模拟了实际中的两个近圆柱状土洞。为便于数值试验，每一粘土层由数层相对独立的组（group）组成，并将模型离散为若干小单元。基岩面大体由两级组成，在土洞下的基岩中用圆柱形空洞模拟了岩溶空洞的存在。几何模型外观及内部土洞见图4-3、图4-4。差分网格的产生由计算机自动生成。
4.4.2.4 力学模型及参数
在建模中所用的力学模型，对于建筑物基础及基岩采用弹性模型（elastic），对于硬塑及可塑粘土层使用摩尔-库仑模型（Mohr-Coulomb）（钱家欢等，1991；雷晓燕，1999）。
在作静力模拟中，FLAC3D所需参数主要包括重力加速度、密度、变形模量、泊松比、内聚力、内摩擦角、抗张强度等。参数来源为郑玉元（1992）的实测资料并参考了《工程地质手册》。计算所需的体积模量及剪切模量，由变形模量及泊松比求得，计算公式如下：
体积弹性模量（bulkmodulus）：

岩溶塌陷机理及其预测与评价研究

剪切模量（shearmodulus）：

岩溶塌陷机理及其预测与评价研究

式中：E——变形模量；μ——泊松比。

图4-3 数值模拟模型外观及差分网格


图4-4 几何模型中的大小土洞

各层参数的选定值如表4-2所示。

表4-2 概化后各层的材料参数表

4.4.2.5 基础与加载
实际中建筑物基础由三个条形基础组成，基础埋深0.8m，用相应的差分单元模拟基础。基础材料为混凝土材料，强度比石灰岩稍次。荷载施加于相应的基础之上，换算压力值为1×105N/m2，模拟建筑物对土层施加的静力。计算中的最大不平衡力收敛良好，平衡比率也达到了计算要求（图4-5）。

图4-5 静力计算中的不平衡力收敛过程

8. 数据处理及建立模型

9.2.2.1 统计量的选取
基于对金刚石/钻石中E型石榴子石包裹体元素含量统计分析来对其产地来源识别，需要预先搜集世界各地已知的前人研究测试的数据，来建立数学模型，以得出产地来源与包裹体元素含量之间的某些联系。表9.1是参与此次统计分析的数据来源及样本数。
参与本文统计分析和绘图等所用的数据，全部来源于该表中对应的文献（附表6）。因此，若下文中无再注明出处或其他特殊说明，其数据均默认来自于该表对应产地的文献，其数字序号也对应相应产地。

表9.1 各产地金刚石石榴子石包裹体电子探针测试数据条数归纳（单位：条）　Table 9.1 Statistics of EPMA test data of garnet inclusions in diamonds from different origins (unit: piece of data)

由于判别分析需要从中筛选出能提供较多信息的变量方能使错判概率变小，因此统计变量的选取尤为重要。包裹体的测试数据包含了数十种元素及其对应氧化物的含量，倘若一一研究，不仅计算量大，计算复杂，而且容易出现重复统计造成较大误差等。在此，作者选定了其中的FeO、MgO、CaO三种组分参数作变量，这样选变量基于如下理由（黄进初，1990）：
（1）Si组分作为石榴子石硅酸盐矿物的主常量组分，不参与此次的统计研究；
（2）Ti、Ni、K、Na、Cr等组分在石榴子石中的含量较低（测试误差大），且测试数据不全（只有某些产地的测试数据，部分产地的测试数据缺失），因此其数值代入统计研究中会引起较大的误差；
（3）FeO、MgO、CaO三个组分是石榴子石中对其种类成分产生主要制约的组分，也是和地幔性质有明显关联性的组分。对归纳的167条E型石榴子石测试数据的预处理显示，各产地FeO、MgO、CaO三个统计变量的数据全，且其组间方差与组内方差比值较大，是各产地间差异性比较大的三种组分参数（其中，Mn2+含量算入Fe2+含量中）。
9.2.2.2 产地归类
对于金刚石/钻石来说，由于其形成环境和条件较为“苛刻”，且世界各产地间由于“历史上”地理位置靠近、幔源性质相近等原因，某些产地间相关包裹体性质具有很大的相似性，仅仅靠一条信息（石榴子石包裹体元素含量统计分析）也许不能区分到具体的每一个产地。为此，本文将先对相似的产地进行归类，研究石榴子石元素含量差异显著的几个代表产地（石榴子石含量差异不显著的产地间将用其他信息来补充区分，本文不详细讨论）。
在此，作者使用主成分综合评价法（陈述云，张崇甫，1995；叶宗裕，2006；阎慈琳，1998），通过将相关统计变量进行主成分分析得到的若干个主成分按线性加权得到一个综合性评价指标，来观察不同产地间的E型石榴子石包裹体地球化学异同。由此，对搜集的各产地E型石榴子石包裹体数据，通过将统计变量FeO、MgO、CaO进行主成分分析，用将所得到的n（1≤n≤3）个主成分按公式9.1提取一个综合主成分：

联合国金伯利进程框架下的钻石原产地研究

其中，Fi为第i（1≤i≤n）个主成分，λi为主成分Fi对应的特征值，λ总为n个主成分的特征值之和。这里取n=2时，其方差累积百分比达99.711%，说明这两个主成分可以很好地综合FeO、MgO、CaO这三个统计变量的信息，且综合主成分值反应产地间的地球化学异同应该具有一定的可靠性。因此，将不同产地综合主成分的平均值作图得到如图9.8。
根据图9.8，将主成分均值相近的产地分为以下四大组，每组内部对应产地E型石榴子石包裹体元素差异较小，而不同组之间差异较明显，可以获得较好的区分度（表9.2）。

图9.8 各产地金刚石E型石榴子石包裹体统计量综合主成分均值图

Figure 9.8 Mean value of comprehensive principle components of garnets inclusions in eclogitic diamonds all over the world

表9.2 产地分组表*　Table 9.2 Groups of diamond origins

表格中的数字序号对应表9.1中的相应产地
9.2.2.3 判别模型的建立
通过判别分析找出各组间的差异性，并建立一个判别模型，作为识别未知产地来源的依据之一。这里使用Fisher判别法，通过坐标变换的方式将数据点投影到另一个坐标系，再用一元方差分析的检验手段将新坐标系中水平差异显著的不同组区分开来，将待判别样本归入离新坐标系中质心最近的组。本文的判别分析过程在统计软件SPSS中进行。
由此，将这4组的FeO、MgO和CaO含量作为统计量，根据表9.2的分组进行判别分析。分析结果部分显示如表9.3所示。
从以上3个表中得到的有用信息如下：
表9.3显示，判别的总判别正确率为67.1%，其中组Ⅰ和组Ⅳ的判别正确率都在80%以上，区分效果较好；但组Ⅱ的正确率仅为50%左右，显示第二组归类样品FeO、MgO和CaO含量的信息与其他几个组之间相关信息的区分度不够明显。
表9.4显示，非标准化的判别方程系数，可以得到一个判别方程组如下：

联合国金伯利进程框架下的钻石原产地研究

其中Ex为判别得分，C为对应物含量。
从表9.5显示，依判别方程，将各组统计量的均值代入可得相应组的质心。若将某个未知来源产地的金刚石E型石榴子石相应FeO、MgO和CaO含量分别代入判别方程9.2组得到的结果E1、E2、E3离哪组的质心距离最近，则认为该金刚石/钻石来源于该产地。图9.9显示，各组质心在同一平面直角坐标系中的位置有显著距离，且各组样本共167条数据作相应转换后的投点归属基本正确，正确率应为67%左右。
图9.10更为直观地显示出不同产地的特征差异：通过四个大组的统计量求氧化物对应的阳离子含量，投Fe-Mg-Ca三元原子百分比图。由于不同组样本数不均，且同组不同产地间仍存在不可避免的部分差异，在以组为单位投点后，为作图的美观性和结果的直观性，每组又再取了一个代表产地的统计量参与对比作图（图9.10）
如图9.10所示，相同产地金刚石/钻石E型石榴子石包裹体Fe、Mg、Ca成分有较好的集聚，而不同产地间又有一定分散的分布，因此具有好的区别性，其中：

表9.3 分组结果　Table 9.3 Regrouping results

*总的判别正确率为67.1%

表9.4 典则判别式函数系数　Table 9.4 Coefficients of Canonical discriminant function


表9.5 组质心处的函数　Table 9.5 Functions at Group Centroids


图9.9 E型石榴子石包裹体产地来源典则判别函数图

Figure 9.9 Canonical discriminant function of garnet inclusion sourcing of eclogitic diamonds

图9.10 金刚石E型石榴子石包裹体Fe-Mg-Ca 原子百分比图

Figure 9.10 Percentage diagram of Fe-Mg-Ca atoms of garnet inclusions in eclogitic diamonds
（1）第一组，加拿大Jericho产地，其平均Mg含量较高，但平均Ca含量较低，Fe含量则分布较散（这里的Fe含量是指Fe2+含量，下同）。
（2）第二组，相应产地的投点则相对分散（图9.9左），这与判别分组结果（表9.2）：组Ⅱ的判别正确率较低相吻合；但其中南非Venetia产地，以其最低的平均Fe含量和最高的平均Ca含量与其他产地有着明显区别（图9.9右）；造成此结果的可能原因将在下文中分析。
（3）第三组，相应产地的投点虽然也有部分分散，但大部分集中在与委内瑞拉Guaniamo产地相近的区域：其特征是各端元含量都居于三个产地之间，这和该组的综合主成分均值也居于所有产地之间结果相吻合。
（4）第四组南非Finsch产地，其特征是平均Fe含量最高，而Ca和Mg含量都相对偏低，与其他组区别明显。
由此可见，不同产地来源的金刚石/钻石E型石榴子石包裹体地球化学性质确存在有较明显的差异性，应该可作为判断未知产地来源的依据之一。