矿床统计预测-实习4a_趋势面分析

实习二、用多元线性回归分析法进行矿床统计预测目的 通过实习，学会使用多元线性回归分析进行矿床统计预测，加深对该方法原理的理解。

要求 （1）根据所提供资料，自己动手完成预测计算的各环节，用控制单元数据建立回归模型，对所提供的各未知单元，评价它们含有小型及以上矿床的可能性；按时提交实习报告。

（2）复习课程“多元线性回归分析”有关内容。

资料 研究区是湖北省某地区一个铁矿成矿带。

为在该区进行矿床统计预测，已将研究区划分为500m ×500m 基本单元408个，并提取了地质变量。

为应用回归分析法，选取一种矿床值（储量乘以某个系数得到的矿床价值）为因变量y 和多个地质因素、找矿标志为自变量(i x )。

本次实习为简便只使用其中两个自变量：1x 为单元磁异常值，2x 为单元中闪长岩体出露面积比。

表2-1（控制单元数据及回归估值表）最左4列列出了实习所用数据。

表2-1控制单元数据及回归估值表方法步骤 第一步：分析研究区内铁矿特征及控矿地质条件和找矿标志，划分基本单元，提取地质变量，地质变量赋值并做适当变换（使因变量尽量服从正态分布、使因变量与自变量之间有线性关系）。

这些工作已经完成（不必重新做），见表2-1的最左4列。

第二步：建立回归方程。

二元回归方程为22110ˆx b x b b y++= （1） 上式中各系数210,,b b b 用最小二乘法确定。

最小二乘法就是使偏差平方和∑=-=ni i i yy Q 12)ˆ( （2） 达到极小而求出210,,b b b 。

上式中n 为观测样品总数。

为此令0,0,0210=∂∂=∂∂=∂∂b Qb Q b Q （3） 将（1）和（2）代入（3），可得到关于210,,b b b 的线性方程组，称正规方程组。

解正规方程组可求出各系数210,,b b b ，从而得到回归方程。

虽然可以采用矩阵形式，但为利于深入细节，本次实习采用分步骤计算的方式。

先计算∑y ，∑1x等值，填满表2-1的最下面一行。

实习1、用证据权法进行找矿远景区预测目的 通过实习，学会使用证据权法进行矿床统计预测，加深对该方法原理的理解。

要求 （1）根据所提供资料，自己动手完成预（2）对计算过程中涉及的计算公式要了解其物理意义；对所涉及各地质变量，要分析了解其地质意义。

（3）复习课程“证据权法”有关内容。

资料 研究区是河北某地区一个北东向复式向斜控制的铁矿集中区。

该区铁矿主要赋存于前铁质来源与火山—沉积作有关，经历了复杂的区域变质（包括混合岩化）和构造变动，矿体多呈大小不等的透镜体状。

方法步骤第一步：分析研究区内控矿地质条件和找矿标志，划分网格单元，提取地质变量（统称为证据层），并将所有地质变量变换为逻辑变量（二值变量），选择控制区（有矿和无矿两类单元）。

在控制单元中统计出各变量存在的单元数（i S ）和含矿单元数（i N ）。

这些工作已经完成（不必重新做），得到表1-1最左边3列。

控制单元总数S =160，其中含矿N =70。

表1-1地质变量（证据层）证据权计算表注：N 表示含有证据层X i 但不含矿的单元数。

第二步：计算各变量的证据权和对比度系数。

证据权分两种，即正权（+i W ）和负权（-i W ）。

它们的计算公式为：)/(/lnN S N N N W i i i -=+)/(1/1lnN S N N N W i i i ---=- (Eq. 1-1)正权和负权分别表示变量与单元含矿和不含矿的关系密切程度。

为表示变量对于单元含矿/不含矿的区分能力，可计算对比度系数（C i ，或称衬度系数），公式为-+-=i i i W W C (Eq. 1-2)根据对比度系数大小可以评价各变量对找矿的重要性。

请根据以上公式，计算填满表1-1，然后填满表1-2。

注意在表1-2中，为节省空间和时间只评价5个变量。

请在每格填写一个变量名（符号）。

表1-2证据层示矿意义评价表第三步：计算各单元的含矿后验概率。

一个变量在任一单元中的证据权为：⎪⎩⎪⎨⎧===-+if ,1if ,i i i i i X W X W W (Eq. 1-3)即若变量在该单元出现，其权为+i W ，否则为-i W 。

矿床统计预测讲义简介矿床统计预测是指通过对已知矿床数据进行统计分析和模型建立，从而对未知矿床进行预测的一种方法。

它是矿床勘探中重要的工具之一，可以帮助矿业公司和勘探者制定科学合理的采矿方案和决策。

本讲义将介绍矿床统计预测的基本原理、主要方法和实际应用，帮助读者了解和掌握该领域的知识和技能。

内容1. 矿床统计预测的基本原理矿床统计预测是基于已知矿床数据的分析和模型建立，通过对已有数据进行统计分析，找出其中的规律和趋势，从而对未知矿床进行预测。

其基本原理包括：•数据收集：收集已知矿床的地质勘探数据，包括地质剖面、岩石样品、地球物理扫描等。

•数据分析：对已有数据进行统计分析，包括数据的中心趋势、离散程度、分布形态等。

•模型建立：根据数据分析结果建立预测模型，包括回归模型、聚类模型、神经网络模型等。

•预测验证：利用已有数据验证模型的准确性和预测能力。

2. 矿床统计预测的主要方法矿床统计预测涉及多种统计学和数学方法，常用的方法包括：2.1. 回归分析回归分析是一种用于探索因变量与一个或多个自变量之间关系的统计方法。

在矿床统计预测中，回归分析可用于确定地质因素对矿床分布的影响程度，并建立预测模型。

2.2. 空间插值空间插值是一种通过已有数据推断未知位置上的值的方法。

在矿床统计预测中，空间插值可用于填补数据缺失的位置，从而得到完整的矿床数据集。

2.3. 聚类分析聚类分析是一种将相似对象归为一类的方法。

在矿床统计预测中，聚类分析可用于将矿床按照地质特征划分为不同的类型，为矿床预测提供参考。

2.4. 神经网络神经网络是一种模拟人脑神经元网络的计算模型。

在矿床统计预测中，神经网络可用于识别矿床数据中的隐藏关系，并建立预测模型。

3. 矿床统计预测的实际应用矿床统计预测在矿业勘探中有着广泛的应用，主要包括以下几个方面：3.1. 矿床评估通过对已有数据的统计分析和模型建立，可以对矿床进行定量评估，包括矿床的储量、品位、开采潜力等指标。

小秦岭金矿田典型矿脉矿化趋势面分析与深部预测小秦岭金矿田是我国西南地区重要的金矿，其矿体主要分布在矿山的东部和北部。

在矿山东部地区，主要是半浅成低温热液多级充填型金矿床；而在北部地区，则主要是深成高温热液单级充填型金矿床。

针对这些矿体，我们可以分析其典型矿脉矿化趋势面，并对其深部进行预测。

针对小秦岭金矿田的半浅成低温热液多级充填型金矿床，其典型矿脉矿化趋势面主要表现为矿脉在深度上有逐渐加深的趋势，同时矿脉宽度也逐渐增大。

此外，矿脉中的石英脉、方解石脉等非金属矿物含量也逐渐增多。

根据这些特征，我们可以初步推断该矿脉的矿化程度逐渐加深，并可能存在与之相伴随的非金属矿物。

对于深部预测，我们可以进一步研究该矿体所处的方位和区域构造特征，以了解其可能的地质条件和地下岩体结构，以便进一步确定开采方案。

而针对小秦岭金矿田的深成高温热液单级充填型金矿床，其典型矿脉矿化趋势面则表现为矿物质种类和含量有明显变化。

矿脉的矿物质种类逐渐由含铜-铜铅矿向含金-金银矿转化，同时含量也逐渐增加。

此外，矿脉的规模也逐渐增大，矿脉内的异常集中度逐渐加强。

基于这些特征，我们可以推断该矿脉由铜铅矿转向含金-金银矿石的可能较大。

对于深部预测，我们需要进一步探究该矿体所处的晶体岩和构造带等地质条件，以确定其形成机理和深部矿化规律，为深部勘探提供科学依据。

总之，了解小秦岭金矿田典型矿脉矿化趋势面对于深入了解矿体的矿化规律和寻找深部矿化提供了指导意义。

针对不同类型的矿体，我们需要采取不同的调研方法和勘探手段，结合地质条件和矿床类型，找到分析方法和策略，实现深部预测的信息全面性和准确性，为金矿勘探成功扩大深入打下了坚实的基础。

MATLAB学习小总结MATLAB英文（Matrix Laborator y）具有强大的数值计算和图形处理功能，以矩阵为基本单元进行运算。

MATLAB在上大学以来就已闻名于耳，又有各种原因却没有机会安装尝试。

在矿床统计预测的实习课中有幸得到老师的讲解，顺利安装后尝试了一些功能，也摸索出了一些知识，现将其做一个小总结如下，当然探索将会继续：Part I 界面认识Command window：执行所有的指令，一次执行一行上的一个或几个语句Workspace：存放变量，同时可以对变量进行保存，重命名，删除等基本操作Current folder:当前工作目录，默认目录为cd('D:\matlab work')（本人安装在D盘）Command history 指令历史窗口，记录执行过的指令和执行时的时间Part II语句形式表达式；变量=表达式（变量的命名规则：字母，数字，下划线，不能以数字开头）Part III基本操作本人学习过程中尝试操作的内容如下类似操作（正在学习）；a=[] b=[1 2 3;4 5 6] c=rand(3) d=eye(3) e=ones(3,3) f=zeros(3) g=diag(c)f=diag(g) a1=c+d b1=a1/2 c1=b1+2 d1=det(a1) f1=inv(a) [v,d]=eig(a1)a2=1:0.1:2 a3=prod([1 2 3 4]) b=magic(3) mean2(b) std2(b) max(b) min(b)size(b) length(b) numel(b) mod(13,6) a=logb x=1:10 y=logx y=log(x) z=exp(x) plot(x,y) z=plot(x,z) z=plot(x,z,'*') plot(x,y,'*') whos save mydata load mydata等等Part IV程序设计含脚本.m文件，程序文件.m,函数文件function现将实习中用到的程序文件、matlab代码运算和制图总结如下：实习一_混合总体筛分Part I 累计概率计算：在command window 中输入matlab code.txt中事先编好的代码，按enter即计算各分组的累计概率，如下：% 全铁品位数据data_TFe = [0.4, 18.6, 38.9, 57.2, 26.7, 51.2, 27.4, 28.5, 28.8, 30.1, 33.9, 39.4, 60.4, 41.1, 47, 47.3, 45.6, 43, 31.5, 34.4, 34.9, 25.4, 50.1, 46.1, 55.8, 45, 57.5, 63.3, 41.7, 56, 37.2, 46.5, 53, 51.9, 48. 4, 40.5, 42.4, 42, 38.5, 38.2, 55.4, 52.8, 27.1, 57.7, 53.5, 56.4, 59. 2, 52.3, 52.8, 50.4, 58.7, 33.7, 59.8, 60.9, 60.1, 51.5, 22.3, 55.3, 59.7, 49.2, 59.4, 54.9, 64.7, 49.3, 35.7, 59.4, 94.2, 56.5, 60.2, 88. 4, 59.8];% 组中值data_mval = 20:4:64;% 数据分组data_xval = [data_mval(1)-4, data_mval, data_mval(end)+4];%频数统计, 忽略特低和特高值count_all = hist(data_TFe, data_xval);count = count_all(2:end-1);% 累积频数count_cum = cumsum(count);% 累积频率freq_cum = count_cum ./ count_cum(end)结果：freq_cum =0.0147 0.0441 0.1176 0.1765 0.2353 0.3529 0.4118 0.5147 0.6618 0.8088 0.9706 1.0000Part II.特异值剔除并绘制直方图运行MF4SPOD_MixDist_Preproc.m文件调用SPOD_mixdist_TFe_Raw中的数据进行特异值剔除，显示剔除特异值后的数据并绘制出直方图如下：MF4SPOD_MixDist_Preproc.m文件：% MF4SPOD_MixDist_PreprocstrPathFileName = 'F:\2学习资料-课件\大三下学期课件\矿床统计预测\1 实习数据与程序\A1_混合总体筛分\SPOD_mixdist_TFe_Raw.txt';data_TFe = load(strPathFileName);% 数据探查figure;plot( data_TFe,'o-');stem( data_TFe,'o');boxplot( data_TFe );hist( data_TFe, 10 );% 剔除特异值v_low = 10;v_high = 80;L = data_TFe > v_low & data_TFe < v_high;data_TFe_Sel = data_TFe(L);close all; % 关闭所有Figure窗口plot( data_TFe_Sel,'o-');% 绘制直方图% [17.5, 20.5, 23.5, 26.5, 29.5, 32.5, 35.5, 38.5, 41.5, 44.5, 47.5, 50.5,% 53.5, 56.5, 59.5, 62.5, 65.5, 68.5]x_hist = 20:4:64;hist( data_TFe_Sel, x_hist );freq = hist( data_TFe_Sel, x_hist );disp( [x_hist;freq]' );结果：20 124 228 532 436 440 844 448 752 1056 1060 1164 220 124 228 532 436 440 844 448 752 1056 1060 1164 220 124 228 532 436 440 844 448 752 1056 1060 1164 220 124 228 532 436 440 844 448 752 1056 1060 1164 2实习三聚类分析导入需进行聚类分析的数据文件：ht_ores_TextID，运行程序文件MPHX_HierarchicalCluster，得到方阵v.pdist 、x.pdist并绘制谱系图程序文件：% MPHX_HierarchicalCluster% X=[% 40, 10, 40, 5, 400, 5, 10, 65, 4, 100, 10;% 40, 15, 40, 5, 85, 250, 10, 15, 8, 30, 10;% 600, 300, 120, 5, 800, 15, 10, 500, 4, 60, 50;% 3000, 650, 120, 350, 400, 5, 10, 500, 18, 30, 50;% 100, 450, 40, 150, 300, 5, 10, 400, 2, 30, 300;% 4000, 4000, 160, 100, 1000, 5, 10, 450, 30, 200, 3000;% 3000, 500, 200, 450, 700, 5, 10, 400, 10, 60, 600;% 60, 400, 40, 250, 200, 5, 10, 150, 1, 100, 100;% 60, 3000, 3000, 60, 1500, 5, 5, 150, 2, 80, 3000;% 200, 3000, 3000, 40, 1500, 5, 40, 150, 2, 60, 1000;% 3000, 150, 60, 40, 300, 5, 40, 250, 2, 100, 10;% 300, 150, 60, 20, 100, 15, 10, 120, 2, 30, 10;% 600, 1000, 100, 10, 100, 25, 10, 70, 6, 100, 3000;% 40, 20, 20, 10, 30, 10, 10, 25, 1, 15, 10% ];mydata = data;mydata = zscore(mydata);ngroups = 3;[ T, X_pdist, Y_linkage, v_cophenet, X4cols_inconsistent ] = ...MFHX_HierarchicalCluster( mydata, 'euclidean', 'complete', ngroups );% 输出个分组的对象for ng = 1:ngroups,objsInGruop = find( T==ng );disp(objsInGruop');end结果：v_pdist =Columns 1 through 154.1342 2.7091 4.5757 2.6123 6.5656 4.3734 1.8761 4.94145.1617 3.5496 1.6400 2.6750 1.9579 4.8160 5.7867Columns 16 through 304.5522 8.10965.9055 4.55256.5987 6.6759 5.4994 3.7505 4.5756 3.8139 3.5372 1.7224 6.1763 3.6672 3.0146Columns 31 through 454.87175.0305 3.7049 2.6236 3.7991 3.2363 3.18005.6308 1.6453 3.85936.1702 6.2835 4.4663 4.1789 4.9514Columns 46 through 604.6701 6.7349 3.2205 2.14095.0627 5.2630 3.8343 1.8785 3.4655 2.4136 5.65666.5909 6.0742 6.9599 6.4353Columns 61 through 757.1931 5.6152 7.6925 3.3149 5.6185 5.9286 4.3258 4.2314 4.6751 4.6975 5.1062 5.3862 3.7241 2.2008 3.0701Columns 76 through 902.56553.5712 6.1976 5.24004.32475.4865 4.8319 5.3364 5.3669 5.5569 3.6759 4.2246 4.0491 2.9721 0.6787Column 913.1998X_pdist =0 4.1342 2.7091 4.5757 2.6123 6.5656 4.37341.8761 4.9414 5.1617 3.5496 1.64002.6750 1.95794.1342 0 4.81605.7867 4.5522 8.1096 5.9055 4.55256.5987 6.6759 5.4994 3.7505 4.5756 3.81392.7091 4.8160 03.5372 1.7224 6.1763 3.66723.01464.87175.0305 3.7049 2.6236 3.7991 3.23634.57575.7867 3.5372 0 3.1800 5.6308 1.6453 3.85936.1702 6.2835 4.4663 4.1789 4.9514 4.67012.6123 4.5522 1.72243.1800 0 6.7349 3.2205 2.1409 5.0627 5.2630 3.8343 1.8785 3.4655 2.41366.5656 8.1096 6.1763 5.6308 6.7349 0 5.6566 6.5909 6.0742 6.9599 6.43537.1931 5.6152 7.69254.37345.9055 3.6672 1.6453 3.2205 5.6566 0 3.3149 5.6185 5.9286 4.3258 4.2314 4.6751 4.69751.8761 4.5525 3.0146 3.85932.1409 6.59093.3149 0 5.1062 5.3862 3.7241 2.2008 3.0701 2.56554.9414 6.5987 4.8717 6.17025.06276.0742 5.61855.1062 0 3.57126.1976 5.2400 4.3247 5.48655.16176.6759 5.0305 6.2835 5.2630 6.9599 5.9286 5.3862 3.5712 0 4.8319 5.3364 5.3669 5.55693.5496 5.4994 3.70494.4663 3.8343 6.4353 4.3258 3.7241 6.1976 4.8319 0 3.6759 4.2246 4.04911.6400 3.75052.6236 4.1789 1.8785 7.1931 4.23142.2008 5.2400 5.33643.6759 0 2.9721 0.67872.6750 4.57563.79914.9514 3.46555.6152 4.67513.07014.32475.3669 4.2246 2.9721 0 3.19981.9579 3.8139 3.2363 4.67012.4136 7.6925 4.69752.5655 5.4865 5.5569 4.0491 0.67873.1998 09 101 2 3 4 5 7 8 11 12 13 146实习四：趋势面分析导入数据文件(以Pb为例)：Juradata.txt% 输入变量x = data(:,1);y = data(:,2);Z = data(:,3);2.调用trend234.m计算函数调用：[B,Z_fit,R2,F] = trend234(x,y,Z,2)Z平均值= 54.630975总离差平方和：SST = 392179.437959回归平方和：SSR = 41547.674570剩余平方和：SSE = 350631.763389解得系数：B = 67.0392 14.5675 -14.6208 -5.0713 6.4953 -2.3474Z_fit =53.3999 65.9024 53.7634 52.7338、、、54.5173 51.7694 52.3233 拟合度R2 =0.1059统计量F =8.36573. 运行trendFig.m绘图Eg:figure1等值线图函数文件：trend234.m,trendFig.m(略)。

MATLAB计算程序1．导入数据文件(以Pb为例)：Juradata.txt% 输入变量x = data(:,1);y = data(:,4);Z = data(:,7);2.调用trend234.m计算函数调用：[B,Z_fit,R2,F] = trend234(x,y,Z,2)Z平均值= 54.630975总离差平方和：SST = 392179.437959回归平方和：SSR = 41547.674570剩余平方和：SSE = 350631.763389解得系数：B = 67.0392 14.5675 -14.6208 -5.0713 6.4953 -2.3474Z_fit =53.3999 65.9024 53.7634 52.7338、、、54.5173 51.7694 52.3233 拟合度R2 =0.1059统计量F =8.36573. 运行trendFig.m绘图Figure1Figure2：Pb原始数据等值线图（二次趋势分析）Figure3.Pb趋势面（二次趋势分析）Figure4.Pb剩余值（二次趋势分析）趋势面分析，是利用数学曲面模拟地理系统要素在空间上的分布及变化趋势的一种数学方法。

它实质上是通过回归分析原理，运用最小二乘法拟合一个二维非线性函数，模拟地理要素在空间上的分布规律，展示地理要素在地域空间上的变化趋势。

通常把实际的地理曲面分解为趋势面和剩余面两部分，前者反映地理要素的宏观分布规律，属于确定性因素作用的结果；而后者则对应于微观局域，是随机因素影响的结果。

分析：本作业以Pb为例进行趋势分析，调用函数作出原始数据等值线图，趋势图和剩余图；根据figure2原始数据等值线图可知该区域Pb高异常主要集中在图中中部及偏南区域，中部一高异常位于森林与草原、耕地的交汇部位，环带明显具有多级异常的特点最高达300ppm；南部凹角亦出现一高异常含量最高达160ppm,另外在西南部亦出现5个环带，形状大小不一，中间最高值80-120ppm。