基于crosta方法的遥感矿物蚀变信息提取操作文档

1、专题四：基于crosta方法的遥感矿物蚀变信息提取1、专题概述蚀变岩石是在热液作用影响下，使矿物成分、化学成分、结构、构造等发生变化的岩石。由于它们经常见于热液矿床的周围，因此被称为蚀变围岩，蚀变围岩是一种重要的找矿标志。遥感地质应用中，近矿围岩蚀变形成的蚀变岩石与其周围的正常岩石在矿物种类、结构、颜色等方面都有差异，这些差异导致了岩石反射光谱特征的差异，并且在某些特定的光谱波段形成了特定蚀变岩石的光谱异常。光谱异常为用遥感图像的异常信息提取提供了理论依据。遥感蚀变异常信息提取的方法有多种，其中主成分分析法提取蚀变信息是相对最为广泛的。主成分分析（PCA：Principal Component

2、 Analysis）是基于信号二阶统计特性的分析方法，由于所获各主成分之间不相关，主成分之间信息没有重复或冗余。多光谱遥感数据通过PCA所获每一主成分常常代表一定的地质意义，且互不重复，即各主成分的地质意义有其独特性。但是由于蚀变矿物形成的影像特征在遥感图像上往往表现得很微弱或不明显，甚至“淹没”在主体色调中。以TM数据为例，通过Crosta方法说明其准则，通过TM1、TM3、TM4、TM5和TM1、TM4、TM5、TM7的波段组合分别进行主成分分析提取铁染蚀变和羟基蚀变信息。由TM1、TM3、TM4、TM5做PCA处理，处理后的某个新的成分可能集中了铁染蚀变信息。对代表铁染蚀变的主成分的判断

3、准则是：TM3的系数应与TM1、TM4的系数相反。由TM1、TM4、TM5、TM7作为输入波段进行主成分分析。对代表羟基和碳酸根离子主成分的判断准则是：TM5系数应与TM7、TM4的系数符号相反，TM1一般与TM5系数符号相同。依有关地物的波谱特征，羟基和碳酸根离子信息包含于符合这判断准则的主成分内。铁染蚀变和羟基蚀变存在于绝大多数成矿岩体中，提取这两种蚀变信息基本可以确定研究区成矿岩石的分布情况。原理及其他方法参考博文：遥感矿物蚀变信息提取方法及ENVI下实现/s/blog_764b1e9d0100so8i.html参考文献：干旱区TM图像蚀变信息

4、提取方法研究、ETM+(TM)蚀变遥感异常提取方法研究与应用-方法选择和技术流程2、处理流程介绍主要包括：图像预处理、PCA分析、选择有效成分、异常切割、蚀变信息处理等步骤（图1）。原始数据经过一系列的图像预处理，包括几何纠正，大气校正，掩膜去背景等；再选择波段进行主成分分析，根据上述的原则进行成分的选择；将选择的成分用异常切割标准处理，突出蚀变信息；并与已知的矿点资料进行叠合分析，验证提取的蚀变结果。在做异常切割时，利用( X + k) 确定异常下限和划分异常强度等级。X是某一成分的统计均值代表区域背景，是该成分的标准差。k取值一般是（13），有了这一标准，切割异常时可以减少主观任意性，并使

5、操作较为规范化。经PCA处理后得到的各成分是互不相关的，理论上可以看成是单一的地物。但是由于大气，传感器、异物同谱等影响，各成分并不只是单一的地物，所以在提取异常时利用阈值切割法。参考ETM+(TM)蚀变遥感异常提取方法研究与应用-方法选择和技术流程TM/ETM+原始数据数据读取辐射定标几何校正大气校正研究区裁剪基准影像蚀变异常等级分布图结果验证1345 PCA分析4铁染异常成分判断数据预处理铁染蚀变信息提取规则裁剪1457 PCA分析羟基异常成分判断异常等级划分异常等级划分羟基蚀变信息提取图 专4-1 流程图1、 详细的处理过程本专题的数据存放在“18-基于crost方法遥感蚀变信息提取”文

6、件夹内。3.1数据预处理第一步：数据读取和定标（1）打开tm原始影像数据选择主菜单- File - Open External File - Landsat-GeoTIFF with Metadata，打开“XXXX LT51450282011249IKR01 L5145028_02820110906_MTL.txt”文件。可以看到ENVI自动进行了波段合成。（2）TM数据辐射定标选择主菜单- Basic Tools - Preprocessing - Calibration Utilities - Landsat Calibration，选择含有多波段的L5145028_0282011090

7、6_MTL.txt，弹出如下对话框，图2，选择保存路径，点击OK。图 专4-2 辐射定标参数设置对话框第二步：FLAASH大气校正（3）辐射量度单位转换选择主菜单- Basic Tools - Band Math ，在对话框中输入转换公式，见下图3 ：图 专4-3 Bandmath 输入单位转换公式 输入转换公式以后，点击Add to List，再点击OK。弹出如下对话框（图4），点击Map Variable to Input File按钮为B1变量赋值，选择经过定标处理后的文件进行单位转换，点击OK。在Variables to Bands Pairings面板中，指定存储路径E:insarc

8、alibration_z.dat，点击OK。转换后的文件自动的添加在波段列表中。具体的操作步骤见下图4：图 专4-4 输入波段组合数据步骤（4）储存顺序调整Flassh大气校正对于波段存储的要求为：BIL，BIP格式，上述计算得到的存储方式为BSQ，在此进行波段存储顺序的转化，具体操作如下：选择主菜单- Basic Tools -Convert Data (BSQ ,BIL ,BIP)；选择上一步处理的结果calibration_z.dat文件，点击OK；弹出如下对话框（图5设置如下，点击OK）：当Convert In Place选Yes时，输出的文件自动的覆盖输入文件。图 专4-5 存放顺序

9、转换（5）Flaash校正参数设置大气校正的前期准备工作完毕，现在进行校正参数的设置：ENVI basic tools preprocessing calibration utilities FLAASH，弹出对话框（图6）：单击Input Radiance Image按钮，选择上一步准备好的辐射亮度值数据。由于经过了单位换算过程，在Radiance Scale Factors对话框中选择Use single scale factor for all bands（Single scale factor：1.000000）。Lat：46.01672145 Log：84.60429840图 专4-

10、6 flaash大气校正参数输入图 专4-7 多光谱设置对对话框图 专4-8 flaash高级设置面板（6）根据上述图中的参数设置（图6,7,8），然后点击ok，运行flaash大气校正。运行结果见图9。注意：可以使用restore按钮导入事先准备好的模板文件template.txt，在“中间结果”文件夹中。图 专4-9 校正结果比较第三步：研究区裁剪一、 工程区裁剪（1） 打开前面做好的文件：reflect.dat，（2） 主菜单-File-Save File As-ENVI Standard，弹出New File Builder面板，（3） 在New File Builder面板中，单击I

11、mport File，弹出的Create New File Input File面板，（4） 在Create New File Input File面板中，选中Select Input File列表中的裁剪数据，单击Spatial Subset按钮，（5） 在Select Spatial Subset面板中，单击Image，弹出Subset by Image对话框（图10），（6） 在Subset by Image对话框中，按住鼠标左键拖动图像中的红色矩形框确定裁剪区域，单击OK，（7） 在Select Spatial Subset面板中，可以看到裁剪区域信息，单击OK，（8） 在Create

12、 New File Input File对话框中，单击OK，（9） 在New File Builder，单击Choose设置输出文件名yanjiuqu.dat及路径，单击OK。完成研究区的裁剪过程。图 专4-10 研究区域裁剪二、 反射率单位转换（说明，为了保证学员练习操作时，采用统一的研究区域，请打开数据文件夹中的“中间结果研究区”文件夹中的yanjiuqu.dat数据。）（1）选择主菜单- Basic Tools - Band Math ，在对话框中输入转换公式b1/10000.0（2）按照“第二步，flaash大气校正中辐射量度单位转换”中操作的步骤，实现反射率的单位转换。B1变量的指定

13、值为yanjiuqu.dat文件。指定保存路径和保存文件名yanjiuqu_z.dat。第四步：掩膜一、消去植被和水体的影响（1）选择主菜单-Vector-Open Vector File，选择river.evf文件（注意：在该练习数据的“中间结果”文件夹中）；（2）选择主菜单-Basic Tools-Masking-Build Mask，选择已经打开的要建立掩膜的文件。（3）在Mask Definition的Options菜单下Import EVFs，选择已经打开的river.evf文件，点击OK，指定输出路径和文件名：mask.dat，点击APPLY，结果文件mask.dat出现在波段列表

14、中。（4）建立反掩膜文件。选择主菜单-Basic Tools-Band Math中输入(b1 eq 0)*1 +(b1 eq 1)*0，指定路径和文件名fan_mask.dat，点击OK（见图11），结果文件fan_mask.dat出现在波段列表中，双击显示（见图12）。 图 专4-11 波段运算中变量赋值图 专4-12 掩膜与反掩膜对比图3.2信息提取一、羟基蚀变信息提取（1）选择主菜单-Transform-Principal Component-Forward PC Rotation-Computer New Statistics and Rotate。（2）在Principal Comp

15、onents Input File中选择yanjiuqu_z.dat文件；点击Spectral Subset按钮，按住ctrl键，选择1457波段，点击OK；在Select Mask Band按钮中选择fan_mask.dat文件（图13），点击OK；再点击OK。（3）在Forward PC Parameters中（图14），选择文件保存的路径和文件名：1457PCA.sta和1457pca.dat，点击OK。图 专4-13 选择掩膜波段图 专4-14 主成分分析参数设置二、确定羟基异常成分选择主菜单-Basic Tools-Statistics-View Statistics File，打开

16、1457PCA.sta文件，弹出图15对话框。根据判别规则，确定第4个成分为含有羟基异常的成分。图 专4-15 统计文件分析窗口3.3异常等级划分（1）选择主菜单-Basic Tools-Statistics-Compute Statistics，选择1457pca.dat文件，在Select Mask Band按钮中选择fan_mask.dat文件，点击OK；再点击OK。（计算统计值，消除水体、植被的影响。图16）（2）在Compute Statistics Parameters中，按默认设置（图16），点击OK。得到PC4的标准差为0.005480。图 专4-16 统计计算窗口（3）等级划分 根据标准差的倍数进行等级划分。一级异常的下限值为3x0.00548=0.01644、二级异常2.5x0.00548=0.0137、三级异常2x0.00548=0.01096。这个K参数的大小根据该区域的已知蚀变带进行确定。l 打开1457pca.dat的第4波段PC Band 4并显示；l 在图像窗口-Overlay-Density Slice对羟基异常结果进行分级（图17）；l 选择File-Output Range to Class Image，可以将分割结果输出为ENVI分类格式。修改编辑数据范围和颜色