化探数据处理成图过程

PAGEPAGE17化探数据处理成图的过程毕武1、2段新力1、2黄显义1、2袁小龙1、2彭仲秋1、2李永华1、21.乌鲁木齐金维图文信息科技有限公司，新疆，乌鲁木齐，8300912.新疆地矿局物化探大队计算中心，新疆，昌吉，8311000前言原测数据水文地质数据处理数据预处理重磁数据处理化探数据处理数学地质数据处理点位数据图打印输出原测数据水文地质数据处理数据预处理重磁数据处理化探数据处理数学地质数据处理点位数据图打印输出数据检查数据整理数据网格化投影变换格式转换电法数据处理数据处理图件绘制输出数据输入图1金维地学信息处理研究应用系统(GeoIPAS)图件制作流程图图像输出文本文件图形文件剖面图等值线图场强图三元图三维立体图点位符号图……1处理步骤化探处理的成果包括：(1)参数统计表；(2)R型聚类分析-谱系图；(3)重复样三层套合方差分析或者重复样合格率计算结果；(4)点位数据图；(5)地球化学图；(6)直方图；(7)组合异常图；(8)综合异常图；(9)远景区划图；(10)单元素异常参数统计(附表册)；(11)异常剖析(附图册)；(12)综合异常登记卡(附表册)。在GeoIPAS系统中，化探数据处理分为以下几个主要步骤：1.1数据检查数值检查，坐标检查，重复样坐标检查。1.2分析处理重复样三层套合方差分析、重复样合格率计算、化探特征参数统计、化探背景值分析、R型聚类分析、因子分析。1.3数据分析数据变换；衬值、累加衬值；数据累加、累乘、比值；异常归一化。1.4网格化离散数据网格化、XYZ数据转网格数据1.5成图点位数据图、彩色等量线图、直方图、组合异常图、单元素异常图、综合异常图、剖析图。1.6单元素异常参数统计1.7综合异常登记卡开始开始数据输入数据预处理网格化成图异常统计数据分析处理结束分析样品与坐标对应数值检查坐标检查重复样提取采样编码统计三层套合方差分析或重复样合格率计算化探特征参数统计R型聚类分析离散数据网格化XYZ数据转网格数据点位数据图直方图地球化学图组合异常图单元素异常图综合异常图剖析图TXTGRD单元素异常参数统计图2化探数据处理成图流程TXTTXTTXTTXTTXT矢量GRD矢量TXT矢量CSV2具体处理过程2.1数据检查我们的数据处理工作从化验室提供的样品分析报告开始，项目要提供坐标和样品对应的分析数据，坐标我们一般取实际工作中的米单位，系统中默认东西向横坐标为X坐标，不加带号，南北向纵坐标为Y坐标，需要时还要提供样品对应的地质编码，我们拿到这个数据后首先进行数据检查，以确保数据中不出现写错、漏填、负数、0、>等字符，如果有这样的情况要找实验室给予纠正。数据准备好后，我们要把数据转换成TXT后缀的文本文件，这就做好了处理前的准备数据工作。2.2分析处理2.2.1重复样三层套合方差分析一般是从分析样品的结果中挑出来重复样的分析值，每组四个样品，按如下顺序排列：11第一次采样第一次分析12第一次采样第二次分析21第二次采样第一次分析22第二次采样第二次分析结果：三层套合方差分析成果--元素:Cu方差来源离差平方和自由度平均离差平方和F-检验值地点间13.224180.73511.152样品间1.252190.0661.482分析间1.689380.044三层套合方差分析成果--元素:Au方差来源离差平方和自由度平均离差平方和F-检验值地点间0.573180.03210.781样品间0.056190.0030.483分析间0.232380.006查表：Fa0.05(18,19)=1.35Fb0.05(19,38)=1.22Cu：11.152>1.351.482>1.22即：F1>FaF2>FbAu：10.781>1.350.483<1.22即：F1>FaF2<Fb说明：一般是F1＞临界值，F2＜临界值。Fa>f0.05()理论值,说明采样与分析误差不显著，此次分析结果可靠，Fb>f0.05()理论值,说明采样与分析误差显著，且以采样误差为主。例如报告中的描述：先列一个表格说明所有元素的F1和F2；重复样三层套合方差分析，39种元素F1值远超过临界值，说明采样与分析误差变化均低于区内地球化学自然变化，数据是可用的。F2值表现为两种情况，Ag、Pb、P、U、MgO、La等6个元素F2值小于临界值，说明采样误差不大于分析误差，以分析误差为主；其余33种元素F2值均大于临界值，说明误差应以采样误差为主。2.2.2地层编码统计从样品相对应的地质图中自动读取地层的颜色号作为编码，实现每个样品的地层自动编码。在GeoIPAS系统中的化探->化探特征参数统计->地层编码统计其处理结果是在输入的数据文件中增加一列，过程是将坐标投影到地质图相应的投影参数上，以相同地层颜色来套取数据，取地层的颜色号为落入此地层的采样点的地层编码，所以要想一次做好，要先对地质图做个预处理，最好是按编码顺序，将颜色号改成相应的顺序号，比如Q的颜色改为1，N的颜色改为2，E的颜色改为3，等等。相对应的出来的结果就是按顺序的地层编码了，可以直接用了，绘制的直方图也是这个顺序了。2.2.3化探特征参数统计化探参数统计成果--元素:Cu地质单元样品数(N)平均值(X)中位数(M)标准离差(S)变化系数(Cv)最大值最小值浓集克拉克值(C)致矿系数(Z)1.算术全区210541.91934.7328.6280.683250.492.491.1038.6162.对数全区21051.5241.5410.3070.2012.3990.3960.042.2553.算术剔除高值全区203939.01234.0222.9730.589108.722.491.02710.7344.对数剔除高值全区21051.5241.5410.3070.2012.3990.3960.042.255注：浓集克拉克值可以选用系统提供默认的，也可以直接输入用户所在不同地区的浓集克拉克值。这个值可以保存为文件，以后导入使用。2.2.4R型聚类分析从谱系图中可以看到，Cr、Ni、Co相关性最好，我们在做组合异常图时可以参考谱系图来做元素的组合，但这个也不能完全作为组合参考，我们通常还是以常用组合方式来做组合异常图，如：Cu-Pb-Zn；W-Sn-Mo；Cr-Ni-Co；Au-Ag-As-Sb。2.3网格化2.3.1工作比例尺与测网密度DZ／T0145-94土壤地球化学测量规范表1工作比例尺与测网密度工作阶段比例尺矩形网格正方形网格点／km2线距(m)×点距(m)点线距(m)区域调查1：2000002000×500100011．1000001000×250～500500～10001～4普查1：50000500×100～250250～5004～201：25000250×50～100125～25016～80详查1：10000100×20～5050～100100～5001：500050×10～2525～501：200020×5～10从上表可以看到，化探在实际工作中一般和物探工作一起工作的常采用矩形网格，但在实际成图中却是用正方形网格的好，不用经过网格化处理，没有数据损失，也符合绘图的正方形网。1:20万化探在全国范围已基本完成，只在新疆等高山区有少量的没有覆盖全，最常采用的是组合样，即：1点／4km2，少量的有单点样的，这主要用转网格数据的方式来将数据做网格化处理。单点样的网格化一般采用500×500到1000×1000的网格距，搜索半径3000×3000到5000×5000。组合样做转网格数据，注意起始坐标为奇数公里，间距为2000，然后要进行一次补空扩边处理，这主要是处理数据成图到图框边。1:5万化探是现在化探普查的主要工作，一般在山区采用自由网的水系样，土壤测量一般用正方形网格采样，1:1万化探作为1:5万化探异常检查的主要工作手段，常采用100×40的矩形网格，少量有100×20、200×40、150×40等加密或者抽稀采样。也有采用物探布设的测网非南北方向的测线。2.3.2下面用几个问题来具体说明。问题1：1:5万土壤测量，网距330×300的怎么设置搜索半径和间距？5万化探，网距500*500的，用GeoIPAS做地球化学图的时候，怎么设置搜索半径?答：a、关于网格距，最好采用正方形网格，网格距为野外采样的点线距的一半，比如野外采样为500×500，网格化的网格距选250×250就可以了，100×40的选为50×40或者50×50。b、网格化的搜索半径选择需要一定的经验，建议你首先看看点位数据图，选择一个既能将所有数据联系起来，又不至于网格后把未采样区填充的参数，一般是网格距的2-3倍。问题2：在进行1:5万化探采样时，乙类异常区采样密度设计为12个点／km2，丙、丁类异常区及其它基岩区采样密度设计为个8个点／km2，网格化时XY坐标的间距、行列数如何选？答：一般1:5万化探点线距250×250米或者500×500米；搜索半径以满足丙、丁类异常区网格后只有个别空白区或者无空区，搜索半径一般1000×1000到2000×2000米；如果对网格化参数设置不太清楚，可以选择“计算网格化参数”；网格化后如果异常比较细碎零乱，可用数据预处理→二维数据滤波处理→周围点平均圆滑处理，X和Y方向点数选择3，原点倍数选择2-4倍。问题3：请教如下问题．1:2.5万土壤测量，线距250m，点距40怎么设置搜索半径和间距？答：搜索半径500*500，间距40*40或者50*40，注意最好与原始点位重叠.问题4：1:1万化探，南北向测线，100×40网格距。答：这样的数据本身符合网格数据的格式，只是以TXT文本格式保存，所以我们可以借助GeoIPAS软件的“数据预处理”→“数据格式转换”→“XYZ数据转为网格数据”功能，这里注意数据坐标的起始值、终止值和网格间距，网格间距就选择100×40，起始值和终止值要保证数据的原始点位都与网格点重合。这里有一个问题，就是有些项目用的是GPS坐标，不一定是与网格点位完全重合，一般我们建议用规则网的坐标。圆滑方法：周围点平均圆滑；圆滑参数：X方向点数：3，Y方向点数：3，原点倍数：3问题5：我的测网不是正南北向是45度方向，用咱们的软件进行网格化后变成正南北的，而且图幅形状不是原来测网的形状。请问这样的问题怎样解决?1:1万化探，北西向测线，100×40网格距。答：这样的数据虽然也是网格数据，但是经过旋转的测线方向与我们GRD所定义的网格间距是有夹角的，最常采用的是离散数据网格化方法。离散数据网格化选择幂指数加权，指数因子：10；搜索范围：圆域，300米，网格化间距：40×40。要用XYZ数据转网格数据的方式，就要知道测线的方位角，采样的点线号，起始坐标；利用点线号转换为实际点线距，然后再把点线距转为网格数据，成图后，统一处理，包括旋转方位角，然后平移到图上的起始坐标位置套图框。2.3.3网格化成图的比较：三角剖分绘制的地化图，有被测线方向拉长的趋势，虽然其能保留原占坐标位的原始值，但成图结果不是很美观实用。距离平方导数加权和幂指数加权网格化，对数据的压制较大，异常线的形态较刻板。Kriging网格化选用了两种搜索范围，一种是圆域：300×300；一种是椭圆域：300×100，并且旋转45度。这两种方法整体上差异不大，只在局部的空区上表现出来不同，对数据的压制较小，异常的形态也保持较好，也是常选用的网格化方法，但其也有一个问题，就是对特高值边上的空区会处理出负值，这个不是我们希望得到的结果。2.3.4经验：在处理1:1万化探数据时，可以根据测线方向和数据坐标点位的不同选择不同的数据处理方法，要遵循的一个原则是尽量使网格化后的坐标点位与原始数据坐标点位接近，网格化搜索的范围越小数据损失越小。对于网格间距以小于等于最小的原始点距为准，如果要保留更多的细节，也可以将网格距减小一倍，如果图面太过零乱，可以用GeoIPAS的二维滤波处理的周围点平均圆滑处理，这里可以加大原点倍数，默认为2，可以选择3-5倍。所有这些处理方法可以先选择一个元素用不同的方法做出结果来比较一下确定效果最好的方法，记住处理参数然后再统一批处理。2.4成图2.4.1图件编制的要求：a、土壤地球化学测量图件分两部分：基础图与推断解释图。b、区域调查和普查工作的图件编制按化探区域调查有关规范执行。c、图件编制必须符合地球化学勘查图式图例及用色标准规定。d、成果报告需做交通测区位置图、实际材料图、等值线图、综合异常图及其它推断解释图。2.4.2彩色等量线图问题1：有些异常为什么不标不上高值点?答：这是因为网格化后成图时，最高的那一个等值线没有套住高值点的坐标位置。对这种情况一般我们只能用手动加注极值点，目前还没有太好的办法解决。程序只搜索那个最高的或者最低的等值线范围内的坐标点，然后从搜索的点中选择最大或者最小的进行标注。加注后是一般删线，保证极值点在最高的或者最低的等值线范围内。有时点刚好在线边上也就不做处理了。加标注时，在MAPGIS工程中，加入点位数据图文件，将其设为关闭状态，先在等值线图上找到没标的地方，然后打开点位数据图，把相应的极值拷贝粘贴过来，然后修改子图+，-。问题2：对于等值线高值的圈圈很密怎么处理?答：对于等值线高值的圈圈很密可以减少或者帛稀特高值的画线间隔，不改变常量范围的等值线。另外，有时为了图面好看也可以处理一下特高值数据，我们在GeoIPAS系统中有这样的处理模块，就是压制一下特高值数据，然后网格化成图，可以视情况选用。使用注意事项：（1）数据横坐标(东向距)不带投影带号。（2）在绘制投影经纬网图框时，在原数据投影参数设置中输入正确的数据坐标单位和投影参数，通常高斯坐标数据的单位为米，比例尺为1。（3）在选择标注原始极值时，必须选择正确的横坐标、纵坐标和元素名称或变量名称，否则极值标注将会出错。2.4.3直方图生成直方图的流程为：首先选择数据文件，程序使用TXT格式数据，其次选定元素和组间隔，当选择绘制全区及各地层单元直方图时需要数据中有编码列，指定编码列，有地层编码文件则装入，没有就手动设置，最后设置剖面绘制起始值等设置，完成直方图绘制。a、数据文件格式为TXT文本格式，数据要包含列标识，数据中必须包含地质编码。b、直方图的最大分组数和起始值可以用户自己选择，并且每个地层编码可以分别进行设置。c、地质编码相当于数字的地质图，将地质体用数字代表，每个样品都对应一个地质时代，当某个时代的样品个数少于30个时可以将其合并到其它地质单元中。d、当用户没有编制地层编码文件时，可临时在编码窗口中挑选编码所在的地层。方法是双击“地层代号”下的空白行。e、使用注意事项：（1）编码列必须选择正确。（2）当修改直方图起始值时，小数点后第二次的正值的尾数为7，负值的尾数为3。f、直方图绘图参数设置（1）选择绘制对应形态曲线：即程序自动根据直方图勾绘是否属于正态分布的圆滑曲线；（2）柱体内填充颜色：即每个柱体内填充相应的颜色；（3）对于X刻度线，用户可以根据实际需要进行设置。2.4.4组合异常图a、我们推荐在做组合异常图时，用单个元素来做，在MAPGIS工程中加入不同单个元素异常线文件来任意组合。b、另外如果要做元素异常参数统计，推荐在MAPGIS中先对单元素异常进行编辑整理压缩保存线文件后用单元素异常参数统计来做。c、生成组合异常图的流程为：首先选择数据文件所在的路径，其次选择元素组合，然后进行异常统计设置，最后绘制图框完成组合异常图。d、组合异常图使用的数据文件为Grd格式的网格数据，并且要求各元素网格化时的参数一致。e、在数据列选择框中选择需组合的元素，每选择一个元素，在异常参数设置表中就会增加该元素，并提供系统缺省的异常下限值和线宽度、颜色等绘制参数，缺省的异常下限计算方法有累频和均值标准差两种方式，两种方式分别提供了弱、中、强3种选择，对累计频率而言弱、中、强分别对应了累计频率的85%、90%和95%，对均值标准差而言弱、中、强分别对应了平均值+2倍的离差、+2.3倍的离差和+2.7倍的离差，如果对系统提供的缺省下限值不满意，可以输入确定的异常下限值；每种元素异常对应的线颜色、线宽、线型均可在异常参数设置表中进行设置。2.4.5单元素异常图单元素异常图使用的数据文件为Grd格式的网格数据，绘制方式有3种：直接用异常下限勾绘异常、采用异常衬度值勾绘异常和采用剩余异常值勾绘异常：a、直接用异常下限值勾绘异常，并以一定的含量等级划分出异常浓度外带、中带和内带。当数据为单一地质单元或数据基本符合正态分布时，按照全区异常下限值的1、2～4、3～8倍划分三个浓度带，勾绘异常外、中、内带。b、采用异常衬度值勾绘异常，当数据为两个以上的多个母体分布时，可依多个母体分别统计出的衬值，统一确定衬度异常下限，勾绘异常。并以衬度异常下限的1、2～4、3～8倍划分三个异常浓度分带。c、采用剩余异常值勾绘异常，以求得普查区区域背景趋势面，然后用原始数据减去该点的背景趋势值所获得的剩余异常值作图。当求得剩余异常值后，以正剩余异常值确定剩余异常下限，并以一定的倍数勾绘正剩余异常及浓度分区。在应用剩余异常值制作异常图时，应对求取地球化学背景趋势面的方法进行试验、研究，以确保圈出的异常具有明确的找矿指示意义。2.4.6剖析图GeoIPAS做异常剖析图有两种方法：地化图切剖析图和GRD数据生剖析图。首先都要准备好切剖析图的矩形框，这个可以用MAPGIS中线编辑绘制矩形线，将异常包括进去。a、地化图切剖析图在GeoIPAS中做，利用地化图切剖析图模块，添加地质图和地化图的MAPGIS文件，然后设置地质图的投影参数，导入裁剪框，设置裁剪框的投影参数，选择文件时，将地质图放在第一个位置，选择结果文件目录，进行排版设置，确定后进行裁剪，裁剪结果就组成一个完整的剖析图MAPGIS工程。推荐在MAPGIS中做，可以用剖析图的裁剪框，做成一个区，然后在MAPGIS地质图编辑中用工程裁剪，把地质图和地化图的裁剪结果文件加入到剖析图工程中，进行排版，组成一个完整的剖析图MAPGIS工程，第一个位置放地质图。b、GRD数据生剖析图推荐在GeoIPAS中用GRD数据生剖析图，地质图可在MAPGIS中用工程裁剪得到，裁剪框也用MAPGIS中线编辑绘制矩形线得到。要注意的是有三处设置投影参数的地方，源数据投影参数就是GRD文件的，生成的结果剖析图的投影参数，以及裁剪框的投影参数，选择排列好剖析的图顺序，可以生成规范要求的三级浓度分带的剖析图。2.5单元素异常参数统计GeoIPAS中有两种方式：a、在组合异常图时，选择上“异常参数统计”输入原始数据文件名以及坐标，选择一个或者多个网格化后的GRD文件，选择异常下限，以及最后一列的统计元素名，然后计算出明码文件的同时还会产生一个结果文件+pSta.csv的属性文件，如结果文件名是Ag(.wal.wat),那么属性文件名AgpSta.csv，这个文件就是和异常对应的参数统计，GeoIPAS同时也在结果线文件上挂接好属性。b、由异常线来做异常参数统计做元素异常参数统计时，我们推荐用组合异常图把每一个元素的异常生成出来，在生图时选择不套图框，不绘比例尺，然后在MAPGIS中编辑对异常线做处理。具体包括：删除小的异常线，闭合被图框或工区边角截断的异常线，连接并闭合被空区截断的线，使元素异常个数与实际统计的物理数一致，然后按编排的异常顺序号修改线的ID号，并压缩保存线文件(MAPGIS编辑子系统中，在左边工程文件栏内，点右键，选择压缩保存工程，就把所有编辑状态的文件都压缩保存了)，然后将这个异常线文件作为输入，在化探->单元素异常参数统计中统计异常参数，注意设置线文件投影参数，这个就是我们生成异常线时的参数，还需要设置数据文件的投影参数。运行后就会有一个CSV的结果，其中包括“元素名-编号”列，如Au-1，就和我们的异常编号对应一致，可以直接当作单元素异常参数统计附表了，同时元素异常线上也挂接好属性了。单元素异常的参数，包括：异常点数、异常面积、异常的平均值、异常的对数平均值（几何平均值）、异常内样品最大值、异常标准离差、异常衬度、异常规模、异常NAP，同时进行了排序。异常编号点数序数异常面积序数异常均值序数几何平均值序数异常极大值序数异常均方差序数异常衬度序数异常规模序数异常NAP序数带数序数排序结果Cu-011453.81983460.5557191.779431969.62196.88729171.513891978.519475.7828261314Cu-022936.40139363.18181.791821892.911213.224151.579518148.384410.1113229Cu-032140.7054951273.41121.8197114105.51045.382121.835251223.5706131.29476132211Cu-04971.93215890.516751.918896195.7244.343832.26292597.605964.37237312Cu-052140.56871652.405251.709522663.442315.6058121.31012257.05472210.745068221322Cu-085101.13022989.12661.921125130936.564862.22815655.523492.518319227Cu-102140.6018811454.09231.72112466.642117.7484111.35225238.4805200.