《计算机地质学》实验指导书

《计算机地质学》实验指导书

《计算机地质学》实验指导书西安科技大学地环系二零零六年八月\n目录实验一预处理与统计1实验二一元线性回归分析7实验三多元线性回归分析11试验四趋势面分析15试验五聚类分析20试验六两类判别分析24试验七贝叶斯多类判别分析28试验八有序地质量最优分割法32\n实验一预处理与统计一、目的：通过完成数据统计和预处理程序的设计和实现及完成算例，掌握统计一组数据的极值、均值、方差、变异系数及进行数据预处理的方法。二、方法概要1、进行统计和预处理的原因、目的和应注意的问题（1）原因原始数据可能有强非对称性，存在孤立值，大多数的统计方法应用原始数据时存在大而且不是偶然的残差等问题，通过改变表达方式，有时可以增强信息的显示，而这种改变不仅需要改变数值的单位，而且可能改变数据的基本测量尺度；（2）目的Ø使变量尽可能为正态分布（如回归分析要求因变量为正态分布，要求自变量和因变量之间具有足够的相关关系）；Ø统一变量的数据尺度；Ø使变量之间的非线性关系转换为线性关系；Ø用新的数目少的相互独立的变量代替相互联系的原始变量；Ø方便用简单自然的方式进行解释；Ø帮助理解数据的特征。（3）注意问题Ø数据范围：只有数据变化范围相对较大，变换才显著；Ø变换是很重要的工作，变换不当则适得其反；所以在认真研究分析的基础上进行，有时要通过多次试验才能找到合适的变换方法；Ø有些行业中，有些强制性变换或习惯使用的变换，工作中应遵循；Ø变换后数据的可解释性也很重要，有时为了不影响解释，宁可不对其转换。三、程序设计框图预处理和统计模块的流程图如图1-1所示。后续的不同多元统计方法对原始数据的分布有不同的要求，但对原始数据都需要进行预处理和统计，所以建议统一输入文件的格式，预处理和统计模块处理后输出文件的格式也按该定义的格式定义。文件的格式建议为：31\n样品个数，变量个数第1个样品号，该样品第1个变量的观测值，第2个变量的观测值，…第2个样品号，该样品第1个变量的观测值，第2个变量的观测值，…第3个样品号，该样品第1个变量的观测值，第2个变量的观测值，……不是所有的原始数据都进行系统的统计和预处理，程序应该提供强大和灵活的人机交互界面，方便用户根据实际情况和需要，选择不同的数据变换方法。图1-1数据预处理和统计模块流程图31\n四、算例（1）请对以下15个样品数据进行数据统计和数据预处理样品变量12345678911.7661153.48137952.625.7015.279.955.829.60115821.7287603.23419252.522.2038.922.401.041.40157831.7154183.68239061.105.2023.352.782.861.36137241.4593923.98246130.3633.5928.213.090.891.58138351.7383053.50571071.904.0014.592.331.041.55128561.7133233.29393453.581.8037.843.231.091.80156471.7517413.48281566.164.0018.594.461.381.61132081.8536373.46121464.265.9014.275.562.273.98123091.6496273.77359255.807.8023.894.193.112.231350101.7209033.39263951.702.6037.193.362.441.001549111.7379083.47131159.203.7030.482.751.180.941479121.6318493.82491856.1210.5017.38.372.074.901233131.7204073.29241656.064.6036.001.171.230.311555141.8287893.49714272.886.0514.493.091.041.281290151.7099483.67559567.905.4015.844.392.962.19128231\n实验二一元线性回归分析一、目的：通过对一元线性回归分析程序设计及完成算例，掌握一元线性回归分析的基本原理和方法。二、方法概要1、原始数据预处理对对原始数据（）(=0，1，…，)数据进行预处理（见实验一）。yz散点图2、数学模型建立（1）作散点图一元线性回归只能处理两个变量之间的关系，它又被称为直线拟合，假设为自变量，为因变量，对预处理后对数据（）(=0，1，…，；其中)如果把各个点数据画在坐标纸上（作散点图），各点近似分布呈一条直线，则可用一元现行回归。（2）选择数学模型（2-1）（3）参数a,b的估计为根据回归方程的得到的因变量的计算值；a，b为回归方程中的系数；x为自变量。由于测定结果中不可避免会有实验误差，因此用最小二乘法的原理估计回归直线中的系数a，b。假设实验得到了n对数据（χi，y）(i=0，1，…，n)并得到了回归方程，则对于x的一系列变量x1，x2，…，xn，根据回归得到方程可得到因变量的一系列计算值每一个实测值yi(i=0，1，…，n)与它相对应的一个计算值(i=0，1，…，n)之间都有偏差，也可称残差，计算公式：（2-2）所有测试数据的残差平方和为31\n（2-3）如果回归方程是合理的a，b为最佳。则所得到的残差平方和应达到最小值。得（2-4）列出回归方程（2-5）3、回归方程显著性检验（1）将自变量n组试验数据依次代入回归方程（式3－8），求出因变量n个回归值（）（2-6）（2）计算，和：（2-7）（2-8）（2-9）（3）计算F统计量（2-10）（4）计算相关系数r（2-11）（5）查表检验给定显著水平α，查表得。若或，则证明回归方程显著；否则，回归方程无实用价值。4、对因变量进行区间估计（1）估计剩余标准值（2-12）31\n（2）求出因变量估计值当自变量取定值时，由回归方程可求得因变量数学期望的估计值。（3）写出预测区间95%置信区间:（）99.7%置信区间：（）三、程序设计框图建议对原始数据以以下格式输入：样品个数，2第1个样品号，该样品第1个变量的观测值，第2个变量的观测值第2个样品号，该样品第1个变量的观测值，第2个变量的观测值第3个样品号，该样品第1个变量的观测值，第2个变量的观测值…首先要对数据进行预处理，处理完成后，形成以上格式新的文件（详见图1－1）。以下步骤如流程图（图2－1）。图2-1一元线性回归分析模块流程图四、算例（一）从某煤矿采集11个煤样，分别测定煤的发热量（）和煤灰分（）含量，获得如下表数据：31\n样品号1234567891011（千卡）8.307.87.77.26.86.25.65.55.004.704.30.030.050.080.100.150.200.270.300.340.400.45试建立煤的发热量（）对煤灰分（）的回归方程，并检验该回归方程的显著性。（二）煤矿脉中13个相邻样本点处某种伴生金属的含量数据如下表：样品号1234567距离x23457810含量y106.42108.20109.58109.50110.00109.93110.49样品号8910111213距离x111415161819含量y110.59110.60110.90110.76111.00111.20试建立y对x的回归方程，并进行显著性检验。31\n实验三多元线性回归分析一、目的：通过对多元线性回归分析程序设计及完成算例，掌握多元线性回归分析的基本原理和方法。二、方法概要设有自变量，因变量y，共做n次实验。若y与间有线行关系，回归方程则为：（3-1）显而易见，只要确定了各回归系数，方程也就确定了。1、建立原始数据矩阵，并进行必要的预处理（1）确定因变量与自变量，形成行列的矩阵：（3-2）（2）进行预处理（详见图1－1）2、数学模型建立设经过预处理的原始数据，表示如下（3-3）由最小二乘法知道，使得全部观察值与回归值的偏差平方和Q达到最小，即使（3-4）根据微积分学中的极值原理，b0，b1，…，bn，应是下列方程组：（3-5）3－5式是求解回归系数的正规方程组，可以写成以下形式31\n（3-6）（3-7）求解正规方程组，得到各回归系数。列出回归方程（3-8）3、回归方程显著性检验（1）将自变量n组试验数据依次代入回归方程（式3－8），求出因变量n个回归值（）（3-9）（2）计算，和：（3-10）（3-11）（3-12）（3）计算F统计量（3-13）（4）计算相关系数r（3-14）（5）查表检验给定显著水平α，查表得。若或，则证明回归方程显著；否则，回归方程无实用价值。4、对因变量进行区间估计（1）估计剩余标准值（3-15）（2）求出因变量估计值当自变量取定值时，由回归方程可求得因变量数学期望的估计值。（3）写出预测区间95%置信区间:（）99.7%置信区间：（）31\n三、程序设计框图建议以实验一的原始数据格式输入，首先要对数据进行预处理，处理完成后，形成实验一的原始数据格式的新文件（详见图1－1）。以下步骤如流程图（图3－1）。图3-1多元线性回归分析模块流程图31\n四、算例（一）、某矿区从18个矿样中测得Cu,Pd,Pt的含量如下所示，试求Pt对Cu和Pd的回归方程，并检验其有无实用价值。编号Cu(%)Pd(g/t)Pt(g/t)编号Cu(%)Pd(g/t)Pt(g/t)10.360.1690.176100.060.0230.04520.050.0350.020110.530.3290.29230.240.2480.262120.060.0650.08040.050.0150.035130.210.2830.23550.220.2160.234140.030.0300.03060.130.0860.092150.330.210.23770.101.3001.669160.010.0300.05080.070.0090.006170.210.2520.2490.160.1600.150180.080.0530.103（二）、由于碳、氢、氧是煤燃烧过程中产生热量的主要元素，对某矿（褐煤、常烟煤、肥煤、焦煤和无烟煤）中取12块煤样，经分析化验后，其发热量（焦耳/克）与Cr(%)、Hr（%）、Or（%）元素的数据如下。问：今后能否不再对该矿煤进行发热量测定，而由元素分析结果对其进行预测？编号Cr(%)Hr（%）Or（%）（J/g）编号Cr(%)Hr（%）Or（%）（J/g）1625.0156.57856.06.08.42706.0207.08883.03.08.53756.5257.29905.05.08.84755.0107.510902.52.58.05785.5157.711923.03.08.56805.24.08.012953.53.58.731\n试验四趋势面分析一、目的：通过趋势面分析程序设计及完成算例，掌握趋势面分析方法原理及工作步骤。二、方法概要1、建立原始数据矩阵，并进行必要的预处理（1）选取应变量选取什么样的变量进行趋势分析，取决于研究对象和研究目的。（2）确定控制点在确定控制点时，一般要考虑以下几个方面：①控制点的代表性、真实性和可靠性。②控制点要具有面性、均匀、随机分布的特点。（3）整理原始数据统计表收集因变量观测值及相应坐标，填入趋势分析原始数据统计表，其格式可参考表4.1。表4.1趋势分析原始数据统计表格式编号控制点坐标变量顺序号原编号横坐标()纵坐标()其中“控制点坐标’’可以是地理坐标，也可以是虚拟坐标系的相对坐标。（4）建立原始数据矩阵假设有n个控制点的地理坐标及地质变量观测值，则原始数据矩阵为：（5）统一量纲为了不改变地质变量的值，仅对控制点的坐标值进行以下变换（先中心化，在均匀化）：其中，，；为常系数，其目的在于使，与具有相同的数量级。例如，在对间隔型（标高、水位等）进行趋势分析时，因为一般有三位整数，故将31\n取作100即可；对比例型数据（煤厚、灰分、显微组分含量等）进行趋势分析时，可视的数量级，为值在中间取值。若自选坐标系，控制点的坐标为相对坐标，可通过适当选取坐标系的单位，达到与地质变量统一量纲的目的，因而不必再施加上述变换。2、趋势方程的参数估计纪正规方程组为：（4-1）解方程组可求出系数，，，…，，从而得到趋势面方程：（4-2）其中，3、趋势与剩余、异常与“噪声”的分离（1）由趋势方程计算各控制点的趋势值，再由观测值、趋势值求得各点的剩余值：（4-3）（2）计算正剩余的平均值（4-4）或正剩余的二倍标准差（4-5）（3）计算异常点和异常值选或作为异常限，记为E。若正异常剩余与异常限之差，则第点为正异常点，为正异常值；若负剩余与异常限之和，则第点为负异常点，为负异常值。4、趋势方程的显著性检验（1）计算趋势平方和、剩余平方和、拟合度及检验统计量（4-6）31\n（4-7）（4-8）（4-9）（4-10）（2）给定置信水平，查临界值表的。若，则认为趋势方程在置信水平下是显著的。5、绘制趋势图和偏差图，并进行地质解释三、程序设计框图要求：针对算例编写，不要求通用性，这样可以少占用一些时间。建议以实验一的原始数据格式输入，首先要对数据进行预处理，处理完成后，形成实验一的原始数据格式的新文件（详见图1－1）。以下步骤如流程图（图4－1）。图4-1趋势分析分析模块流程图31\n四、算例（一）、在15*15km2范围内，选择15个观测点测量某沉积物的粒径（坐标单位：km；粒径单位：mm）编号坐标粒径（）mm编号坐标粒径（）mm横坐标()km纵坐标()km横坐标()km纵坐标()km1231.998111.322102.3101081.232131.11111131.44412.6121231.75582.2131261.865141.81412101.27733.51515131.08763.11618161.4要求进行趋势面分析，具体如下：（1）建立合适的面趋势面方程；（2）求出各点的；（3）计算F值并进行检验；（4）绘趋势图和偏差图；（5）进行地质解释。（二）、对下表中煤层底板标高数据（）进行趋势面分析。要求：（1）建立合适的面趋势面方程；（2）求出各点的；（3）计算F值并进行检验；（4）绘一次趋势图和一次偏差图，比例尺为1：5000；（5）进行地质解释。编号钻孔号1.`11a`247.5987.5-38.0-88.0-38.0-38.02.`12a`217.5886.0-60.0-110.0-60.00-60.03.`13`155.0611.0-105.0-155.0-105.0-105.04.`21`540.01125.0-13.5-63.5-13.5-13.55.`22`486.5902.5-54.5-104.5-54.5-54.56.`23`431.5673.0-95.5-145.5-95.50-95.531\n1.`24`367.5424.0-141.5-191.5-141.5-141.52.`25`282.587.5-204.0-245.0-204.0-204.03.`31`，777.51049.5-26.5-76.5-76.5-76.54.`32`，717.080.0-71.0-121.0-121.0-121.05.`33`，656.5536.5-119.0-169.0-169.0-169.06.`34`，591.0284.0-166.0-216.0-216.0-216.07.`35`，559.067.5-207.0-257.0-257.0-257.08.`41`，1055.01121.5-63.5-13.5-63.5-63.59.`42`，986.0849.0-113.5-63.50-113.5-113.510.`43`，922.0576.0-160.0-110.0-160.0-160.011.`44`，858.0341.0-203.5-153.5-203.5-203.512.`45`，804.5119.0-254.0-195.0-245.0-245.013.`51`，1276.0975.0-90.0-40.0-140.0-40.014.`52`，1214.5722.0-134.0-84.0-184.0-84.015.`53`，1154.5490.0-176.0-126.0-226.0-126.016.`54`，1060.0128.0-240.5-190.5-290.5-190.517.`61`，1575.01107.5-67.0-17.0-117.0-17.018.`62`，1502.5826.0-116.0-66.00-166.0-66.0019.`63`，1456.5624.5-151.0-101.0-201.0-101.020.`64`，1385.0388.0-196.0-146.0-246.0-146.021.`65`，1329.5171.0-234.9-184.9-284.9-184.931\n试验五聚类分析一、目的：通过对聚类分析简单程序设计及算例计算，掌握聚类分析方法步骤及其数学原理。二、方法概要1、合理地选择变量和样品，构造原始数据矩阵设有n个样品，m个变量，则原始数据矩阵其中，为第i个样品第j个变量的观测值（i=1,2,…,n;j=1,2,…,m）。2、数据预测理（统一量纲变换）通常，对服从正态分布的变量进行标准化变换，对服从对数正态分布的变量进行对数化变换，而变量不服从正态分布，或各变量服从不同的分布时，则采用正规化处理（详细处理过程见实验一）。3、计算聚类统计量一般在进行R型聚类时，计算相关系数；进行Q统计量型聚类时，计算相似系数或距离函数，其中以斜交距离函数的分类效果较好。事实上，各种方法各有所长，在计算中可选用多个统计量分别进行计算，然后结合实际情况，选择合适的方法。（1）R型聚类采用相关系数或相似系数，按下式计算相关系数：（5－1）其中，，分别为第j个变量和第k个变量的样本均值，分别为第j个变量和第k个变量的样本方差。当原始数据经过标准化数据后，由于，，（6－1）式可简化为（5－2）（2）相似系数（S）31\n（5－3）（3）欧氏距离系数（D）（5－4）（4）斜交距离系数（）（5－5）：为第j个变量与第l个变量的相关系数。（5）误差平方和增量（ΔE）若P点群与q点群合并为t点群，误差平方和增量为Δ,则（5－6）4、划分类别（点群归并）可用一次计算或逐步计算聚类法对变量或样品进行分类。两种方法比较而言，后者较前者较为合理一些，但计算量比前者约大一倍。（1）对相关系数（R）和相似系数（S）,用逐步计算法或加权平均迭次推算法①修改变量数据Q型：i，k合并，则为（5－7）式中，分别为第。R型：j,k合并，则（5－8）式中，分别为第j,k点群中变量的个数。②加权平均迭次推算法若P点群与q点群合并为t点群，则t点群与其他r点的群的相似性系数为：31\n（5－9）（2）对于D、和ΔE采用刷新方程迭次推算法①D、的刷新方程（含：）若P,q合并为t，则t与其他r点群的为：（5－10）式中：分别为p,q点群中的样品个数。②误差平方和增量（ΔE）的刷新方程。若p,q点群合并为t点群，则t与r点群合并的ΔE为（5－11）式中，分别为t,r,p,q点群中样品的个数。5、整理联结表，绘制谱系图6、对聚类结果进行解释三、程序设计框图要求：对算例设计程序，不要求程序具有通用性。聚类分析程序框图（见图5－1）31\n图5-1聚类分析模块流程图四、算例从某地超基性岩石的某些样品中，得到与矿化度有关的一些元素的光谱分析数据。依次为：Ni，Co，Cu，Cr，S，As的原始数据，进行聚类分析。NiCoCuCrSAs190323287442782177510642737926273944416061150136160317584124003140209381635864258234541044143371431\n试验六两类判别分析一、目的：通过对两类判别分析程序设计及完成算例，掌握两类判别分析基本原理和方法。二、方法概要1、整理原始数据，构成数据矩阵，并进行必要的预处理设有n个样品，每个样品测得p个变量的数值，且在这n个样品中，已知有个样品属于A类，个样品属于B类，个样品的归属待定，但必然属于A类或B类。此时，原始数据矩阵如下：（6－1）为消除量纲的影响，对矩阵进行统一量纲变换，变换后的矩阵仍以表示。2、利用A类和B类样品的变量观测数据（经过统一量纲变换），建立判别函数（1）计算组内平均值和组内方差；（6－1）（6－3）（2）计算各变量的值(j=1,2,…,r)（6－4）（3）挑选变量根据Ij的大小，从r个变量中，挑选出P个变量（如）参加建立判别函数。3、建立判别方程（1）计算P个变量的组内方差和协方差31\n(i,j=1,2,…,P)（6－5）（2）解正规方程组（6－6）求出C1，C2…..Cp,获得判别函数4、计算分界值（），得出判别规则（1）计算判别值、类平均值和总平均值（k=1,2,…,nB）（6－7）（6－8）（2）设：（6－9）（3）总结出判别规则若有，则判别规则如下：若待判样品的判别值，则将其归入A类；若待判样品的判别值，则将其归入B类；5、待判样品判别归类计算未知样品得判别值，其中(j=1,2,…,P)若待判样品的判别值，则将其归入A类；若待判样品的判别值，则将其归入B类；6、显著性检验设：（1）计算F值31\n（6－10）在给定下，查F分布表得，如果，则判别函数有效，否则，判别函数无效。（2）计算回代正确判别率（）者为正确，共有个，者为正确，共有个,正确判别率：（6－11）三、程序设计框图图6-1二类判别模块流程图31\n四、算例今获得有关内陆泥炭和滨海泥炭得锶（Sr）和钡（Ba）的化验数据如下：已知类样品号锶（Sr）钡（Ba）Sr/Ba滨海泥炭10.00120.00011220.00300.0005630.00030.00021.540.00520.00182.8850.00020.00021内陆泥炭10.00090.00220.4120.00070.00130.5430.00250.01100.2340.0280.00600.47未知个体10.0010.00052.0020.0050.000657.6930.00070.00023.50试用两类判别分析判断（未知个体）煤层的成因类型。31\n试验七贝叶斯多类判别分析一、目的：通过对多类线性判别分析程序设计和完成算例，掌握贝叶斯多类判别模型原理和工作方法步骤。二、方法概要1、原始数据获取设有G类母体，从每个母体中取得个样品，每个样品测得个变量，则原始数据为：（；；）（总样品个数）2、准备工作（1）计算诸变量的类平均值和总平均值（7－1）（其中；；）或（7－2）（2）计算组内离差矩阵（）和总离差矩阵（T）（7－3）（7－4）()（7－5）（3）求W，T矩阵的逆矩阵（W-1,T-1）及行列式值3、判别分类（1）计算判别系数和及先验概率（7－6）(；)31\n()（7－7）,()（7－8）（2）检验P个变量的判别效果用威尔克斯准则来检验，即检验p个变量对于区分G个母体的能力。计算（7－9）其中：当时当时在给定下，查F分布表得，如果，则判别函数有效，否则，判别函数无效。（3）计算未知个体的判别值：（）（7－10）（4）对未知个体类别分类若()，则将X样品划归第个母体。4、计算后验概率（7－11）式中：31\n5、正确判别率估计设对已知类型n个样品判别归类后，有m个样品归类正确，则正确判别率（7－12）三、程序设计框图图6-1贝叶斯多类判别模块流程图31\n四、算例（一）、某煤矿矿井开采A、B、C三个煤层，由于断层破坏，掘进巷道遇到一层没不知属于哪一层，从而影响掘进工作正常进行。试用判别分析解决该煤层对比问题。为此，取若干煤样，经过化验获得如下数据：煤层样品SAl2O5C2O已知类型A17.8323.352.7427.5823.203.1538.5123.894.1948.3124.004.32B14.7338.922.4025.1237.843.2335.7837.193.3646.1737.403.30C16.5614.592.3327.7815.844.3937.2814.959.3247.3213.943.33未知个体X14.5639.021.36试建立多类线性判别函数，并对未知个体进行判别，求出X属于各母体的后验概率。31\n试验八有序地质量最优分割法一、目的：通过对有序地质量最优分割法简单程序设计及完成算例，掌握该法的数字原理与方法。二、方法概要1、取得原始数据设有N各有序样品，每个样品测得P各变量，则有原始数据矩阵（8－1）2、数据正规化（8－2）3、计算段直径矩阵D（8－3）其中：得（8－4）4、由D矩阵计算全部分段的组内离差平方和，并求出各段的最优分割。（1）最优二段分割1）由D矩阵对每个m=N,N-1，……，2计算相应的组内离差平方和。(j=1,2,…,m-1)（8－5）31\n1）找出最优二段分割点（）（8－6）2）N个样品的最优二段分割为：{x1,x2,…..xa1(N)}{xa1(N)+1…..xN}（8－7）最优二段分割点为a1(N)。（2）最优三段分割1）由D和最优二段分割结果，对每个m=N,N-1,…..3计算相应的组内离差平方和。（j=2,3,…,m-1）（8－8）2）确定最优三段分割及分割点　（）（8－9）N个点的最优三段分割为：（8－10）最优三段分割点为：（3）最优K段分割1）由D和K－1段分割计算的结果，对每个m=N,N-1,….k分别计算相应的组内离差平方和。(j=k-1,k,…,m-1）（8－11）2）确定最优K段分割及分割点（）（8－12）N个样品的最优K段分割为：（8－13）最优K段分割点为：5、给定K值或绘制W-K曲线，确定最终分段段数。31\n三、程序设计框图四、算例(一)、某煤矿所采煤层牌号为主焦煤。在煤巷中见一火成岩墙侵入煤层，致使煤质发生变化，为弄清煤质变化情况，从火成岩附近，每隔0.5m依次取一煤样，现有6个有序煤样的镜煤最大反射率测定数据为试将这些有序数据进行最优三段分割，并说明岩墙对煤质有何影响？31