矿床统计预测-实习2A-证据权法

1、实习2 用证据权法进行找矿远景区预测目的通过实习，学会使用证据权法进行矿床统计预测，加深对该方法原理的理解。要求（1）根据所提供资料，自己动手完成预测计算的各个环节，按时提交实习报告。（2）对计算过程中涉及的计算公式要了解其物理意义；对所涉及各地质变量，要分析了解其地质意义。（3）复习课程“证据权法”有关内容。层次A 手工计算、MS EXCEL、MATLAB辅助计算B 使用GeoDAS矿床资源定量预测GIS软件A 手工计算、MS EXCEL、MATLAB辅助计算资料 研究区是河北某地区一个北东向复式向斜控制的铁矿集中区。该区铁矿主要赋存于前寒武纪变质岩系中，铁质来源与火山沉积作有关，经历了复杂

2、的区域变质（包括混合岩化）和构造变动，矿体多呈大小不等的透镜体状。方法步骤第一步：分析研究区内控矿地质条件和找矿标志，划分网格单元，提取地质变量（统称为证据层），并将所有地质变量变换为逻辑变量（二值变量），选择控制区（有矿和无矿两类单元）。在控制单元中统计出各变量存在的单元数和其中的含矿单元数。假定这些工作已经完成（不必重新做），得到表1-1最左边3列。数学记号说明：S: 控制区，控制区S单元总数为N(S) = 160D: 控制区内矿床(点)，控制区内D所占单元数N(D) =70Xi: 第i个证据层，控制区内Xi所占单元个数N(Xi)Di: Xi包围的矿床(点)，控制区内Xi包围的矿床点(D)

3、所占的单元个数N(Di)表1-1地质变量（证据层）证据权计算表证据层（Xi）N(Di)N(Xi)N(Xi)-N(Di)Wi+Wi-CiCi序次X1 黑云母混合岩片麻岩类64128640.251X2 含紫苏辉石混合岩片麻岩2143220.205X3 含石榴子石混合岩片麻岩3653X4 基性岩残留体1422X5 麻粒岩残留体5186X6 片理平均倾角501032X7 片理平均倾角50-601956X8 片理平均倾角60-702954X9 片理平均倾角701219X10 北东向断裂61102X11 东西向断裂5094X12 北西向断裂2237X13 南北向断裂3456X14 磁异常对数和10327X

4、15 磁异常对数和10-301139X16 磁异常对数和30-501021X17 磁异常对数和50-1002241X18 磁异常对数和1002433X19 磁异常均方差0.5135X20 磁异常均方差0.5-12268X21 磁异常均方差14758第二步：计算各变量的证据权和对比度系数。证据权分两种，即正权（Wi+）和负权（Wi-）。根据，有计算公式为： (Eq. 1-1)正权和负权分别表示变量与单元含矿和不含矿的关系密切程度。为表示变量对于单元含矿/不含矿的区分能力，可计算对比度系数（Ci，或称衬度系数），公式为 (Eq. 1-2)根据对比度系数大小可以评价各变量对找矿的重要性。请根据以上公

5、式，计算填满表1-1，然后填满表1-2。注意在表1-2中，为节省空间和时间只评价46个变量。请在每格填写一个变量名（包括地质说明）。表1-2证据层示矿意义评价表从大到小排序对找矿最有指示意义的变量（根据C判断）出现时对找矿最有利的变量（根据W+判断）不出现时对找矿最有利的变量（根据W-判断）12345从大到小排序对找矿最有指示意义的变量（根据C判断）出现时对找矿最有利的变量（根据W+判断）不出现时对找矿最有利的变量（根据W-判断）1X3 含石榴子石混合岩片麻岩X3 含石榴子石混合岩片麻岩X5 麻粒岩残留体2X9 片理平均倾角70X9 片理平均倾角70X8 片理平均倾角60-703X10 北东向

6、断裂X13 南北向断裂X10 北东向断裂4X18 磁异常对数和100X18 磁异常对数和100X18 磁异常对数和1005X21 磁异常均方差1X21 磁异常均方差1X21 磁异常均方差1第三步：计算各单元的含矿后验概率。一个变量在任一单元中的证据权为： (Eq. 1-3)即若变量在该单元出现，其权为，否则为。对每个单元，所有变量证据权总和F为（省略了W0）： (Eq. 1-4)式中p为变量数，单元的后验概率P(D|X)为： (Eq. 1-5)根据该后验概率大小可评价该单元的找矿有利性。限于时间，本次实习中只计算部分单元（16个）的后验概率。请填写表1-3（单元数据及后验概率计算表）。第四步：

7、预测评价根据单元后验概率大小，评价各单元找矿有利性；根据控制单元的后验概率，选择合适的后验概率临界值，筛选找矿有利单元，圈定找矿远景区；分析有利单元的地质情况，参考后验概率大小，按照成矿预测工作的一般要求，对远景区进行分级。限于时间，本次实习略去详细工作，请根据表1-3，做出合适的结论，填写到该表后面所留空白处。由表1-3，对结果进行分析可知： 由此确定指示找矿前景的后验概率为： 认为找矿有利的未知单元为： 表1-3单元数据及后验概率计算表单元含矿X1 X2 X3X4 X5 X6X7 X8 X9 X10X11X12X13X14X15X16X17X18X19X20X21FP(D|X)111011

8、001000110100000012011010001100100010100311010000101110001000104010101000101000010001500110101000110010000106010010100001010001001070100100101010100001008110101000101010010000191101000101010000011001000111000100110010001001110011010010100010001012011000001101010000010131110010000011001001001401011010

9、00101001000011500110001010100001001016010100010011001000001注：含矿情况，1表示已知有矿，0已知无矿，空白为未知。各变量值1和0分别表示出现和不出现。选用最佳的证据层组合进行计算。B 附录B1. 使用MS EXCEL本次实习因涉及数据及计算量较小，可使用Excel。（1） 表1-1中W+、W-和C的计算：将该表拷贝到一EXCEL工作表，输入公式计算，再将结果拷回到表1-1。例如计算Wi+,公式为形如：=ln(X/70)/(Y/90)，在单元格中键入=号进入公式计算模式。（2） 表1-1中对变量按W+或W-或C排序：将变量一列及排序标志（

10、如C）一列拷贝到一个EXCEL工作表，然后执行排序命令，再将排好的序号序列拷贝到表1-1。（3）表1-3中F及P(D|X)的计算: Step 1: 计算每个单元每个变量的W值：首先将表1-3拷贝到一个工作表，称数据块A；再将表1-1中W+、W-两列拷贝到数据块A右边，称数据块B，注意相应行对齐；再将数据块A拷贝一份到B的右边，称数据块C。Step2: 在数据块C的第一行第一列单元格中输入公式 =IF(Q, X, Y)。其中Q为数据块A的第一行第一列单元格名称，表示逻辑条件真(=1)或假(=0)，与各但与中证据层的取值正好是吻合的；Y和Z分别为数据块B的第一行两列单元格名称，即W+和W-的取值，

11、并都用$固定列号；执行该公式，Q为真，将得到X值，Q为假，则得到Y值，于是，数据块C第一行第一列的值成为相应的W+或W-值。将该第一行第一列单元格公式类推应用，向下，向右拖动，算出数据块C中所有单元格的W。Step 3: 计算F和P(D|X)：在数据块C最下面一行的下面，第一列单元格中输入公式=SUM(X)，这里X表示数据块C第一列的单元格范围。执行该公式得到第一个F。向下拖动复制，算出所有F。用类似方法算出所有P(D|X)。【注意：为叙述方便，这里计算的F值没有考虑W0，请自行加上】Step 4: 最后将结果拷贝到表1-3相应的行中。附录2. 使用MATLAB【MATLAB脚本m文件：wof

12、e_ex.m】% 证据权计算实习定制脚本(1)启动MATLAB，新建m文件; (2)复制本代码段拷贝到m文件中, 保存为文件wofe_ex.m; (3)填写正确的4处 % TODO: 提示语句, 完成计算或显示功能; (4)Run(F5)执行之, 保证代码正常运行及获得预期正确结果.% created by 2013-11-01% % S: 控制区，控制区S单元总数为N(S) = 160% D: 控制区内矿床（点），控制区内D所占单元数N(D) =70% Xi: 第i个证据层，控制区内Xi所占单元个数N(Xi)% Di: 第i个证据层，控制区内Xi包围的

13、矿床点(D)所占的单元个数N(Di)% 证据层（Xi）% X1 黑云母混合岩片麻岩类% X2 含紫苏辉石混合岩片麻岩% X3 含石榴子石混合岩片麻岩% X4 基性岩残留体% X5 麻粒岩残留体% X6 片理平均倾角50% X7 片理平均倾角50-60% X8 片理平均倾角60-70% X9 片理平均倾角70% X10 北东向断裂% X11 东西向断裂% X12 北西向断裂% X13 南北向断裂% X14 磁异常对数和10% X15 磁异常对数和10-30% X16 磁异常对数和30-50% X17 磁异常对数和50-100% X18 磁异常对数和100% X19 磁异常均方差0.5% X20

14、磁异常均方差0.5-1% X21 磁异常均方差1 % 已知模型单元NS = 160;ND = 70;NDi = 64, 21, 36, 14, 51, 10, 19, 29, 12, 61, 50, 22, 34, 3, 11, 10, 22, 24, 1, 22, 47; % 各列表示 X1, X2, ., X21NXi = 128, 43, 53, 22, 86, 32, 56, 54, 19, 102, 94, 37, 56, 27, 39, 21, 41, 33, 35, 68, 58;Evidences_Class = 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3,

15、 3, 3, 4, 4, 4, 4, 4, 5, 5, 5; % 证据分类，本实习中共有5类，分别用1,2,3,4,5来标记。 % 待求预测单元，行：单元，列：证据出现情况。 各列表示 X1, X2, ., X21Units_Pred = .1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1;0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0;1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0

16、, 1, 0;0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1;0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0;1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0;0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0;1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0,

17、1, 0, 0, 0, 0, 1;1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0;0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0;1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0;0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0;1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1,

18、 1, 0, 0, 1, 0, 0, 1, 0, 0;0, 1, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1;0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0;0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1; % 证据权模型计算N_Evidences = numel( NDi );L = N_Evidences = numel( NXi ) & N_Evidences = num

19、el( Evidences_Class );if( L = 1 ), disp(Error Size); return; end minus_NXiDi = NXi - NDi;minus_NSND = NS - ND;r_plus_NDi_ND = NDi ./ ND;r_plus_NXiDi_NSND = minus_NXiDi ./ minus_NSND;W_plus = log(r_plus_NDi_ND ./ r_plus_NXiDi_NSND); % W+W_minus = log(（1-r_plus_NDi_ND .)/（1- r_plus_NXiDi_NSND）);% TODO

20、: W-C_values = W_plus - W_minus;O_D = ND/(NS-ND) ; % TODO: O(D)W0 = log(O_D); % 显示模型单元的计算结果Rlt_Table1 = 1:N_Evidences; Evidences_Class; NDi; NXi; minus_NXiDi; W_plus; W_minus; C_values; % 表格1计算数据idx_W_plus = 6; idx_W_minus = 7; idx_C_values = 8;Rlt_Table1_Dscending,idx_Descending = sortrows(Rlt_Tabl

21、e1,2 -idx_C_values); % 按证据层类别和对比度C值进行排序，-号表示该列降序gscatter(Rlt_Table1(:,1), Rlt_Table1(:,end), Evidences_Class); xlim(0,N_Evidences+1); grid on % 绘制分组散点图，填写Ci分组排序值disp(Rlt_Table1 - 显示表格1计算数据 (第2列为证据层分类，最后三列为 W+ W- C); disp(Rlt_Table1);disp(Rlt_Table1 - 按证据分类排序 (第2列为证据层分类，最后三列为 W+ W- C); disp(Rlt_Table

22、1_Dscending); % 自动提取各个分组中对比度C值最大的证据层作为有利证据层，L_idx_max_C = diff( 0 Evidences_Class ) = 1; % 获得分组降序值最大值所在位置W_Valid = Rlt_Table1_Dscending(L_idx_max_C,:); % 获得证据层序号以及W+,W-等信息disp(W_Valid - 有利证据层获取 (第2列为证据层分类，最后三列为 W+ W- C); disp(W_Valid); % 计算预测单元F值F_post和后验概率P_postUnits_Pred_Valid = Units_Pred( :, W_Valid(:,1) );% 获得预测单元中的有利证据层信息nUnits, nP = size(Units_Pred_Valid);F1_P1_P0_post = zeros(nUnits,3); % F P1(W0不计入)