1∶50 000水系沉积物地球化学勘查采样点位自动设计的方法研究

廖国忠， 张 伟， 吴文贤， 梁生贤

(中国地质调查局 成都地质调查中心，成都 610081)

0 引言

地球化学勘查是一种科学的找矿方法，在找矿实践中已经取得了许多重大成果。在区域地球化学扫面发现的成矿预测区内，基于1∶50 000水系沉积物测量圈定找矿靶区,已成为地质矿产普查的常用手段，且在我国大部分地区取得了较好的应用效果，尤其在金矿勘查方面已成为首选手段[1-3]。自该方法引入国内以来，在众多找矿方法中，一直处于举足轻重的地位，得到广泛地应用。

水系沉积物测量是以水系沉积物为采样对象所进行的地球化学勘查工作,采样点位的设计是否合理，是否能最大限度控制工作区汇水面积，直接影响着水系沉积物地球化学勘查的效果。

当前，水系沉积物地球化学勘查采样点位的设计，主要是根据地形图上的水系和等高线信息，手动勾绘水系分布图，然后再根据水系分布图手工确定采样点位。长期以来都是纯手工作业的方式，工作强度大，对工作人员工作经验要求较高，即使有相当工作经验的不同工作人员所设计的采样点位也有较大差别。让工作人员从繁重的纯手工设计采样点位的工作中解放出来,是笔者力图解决的问题。

前人针对这一纯手工作业方式进行了些许改进：王海鹏等[4]提出的“MAPGIS条件下完美实现不规则网化探采样点自动编号”； 陈可忠等[5]提出的“化探采样点自动编号及航点自动生成的方法”；冯伟华等[6]提出的“基于VBA实现设计点位编号和坐标提取”。通过研究这些文献发现，目前研究仅局限于对采样点点位进行自动编号，尽管此类研究对提高工作效率有了一定的改观，但是对于采样点位如何自动合理的布置这个关键问题完全没有涉足。

笔者提出了基于GIS自动提取水系的1∶50 000 水系沉积物采样点位自动布置的工作方法，极大地提高了工作效率，具有重要的实用价值。

1 研究思路

通过分析研究水系沉积物采样点定位的工作流程可知，最基本的步骤是通过已知水系和地形地貌的特征，勾绘出水系分布图。因此水系地提取是该项工作最基本的步骤。

水系自动提取是地理信息系统GIS和水文科学相结合的交叉研究领域，其中水文特征提取是一个较为常见的功能，已经在生产中得到了广泛地应用[7-9]。因此要实现水系沉积物地球化学勘查中最基础的水系自动提取，引入水文特征提取这一功能便可解决。

水系提取成功后，继续提取水系中的特征点(水系入口，出口，中间点等)作为备选采样点。根据《地球化学普查规范》，1∶50 000水系沉积物地球化学勘查工作中，采样点布设密度为4个点/km2～8个点/km2；而通过提取所获得的备选采样点的密度有可能远远大于4个点/km2～8个点/km2。因此，必须根据《地球化学普查规范》建立一个评价体系，选取最佳的备选采样点作为最终的采样设计点位。将筛选的采样设计点位根据《地球化学普查规范》要求自动编号。

2 采样点位自动设计方法的实现

采样点位自动设计方法的实现需要以下几个步骤:

2.1 已知水系限制下的水系提取

以获取的1∶50 000数字地形图为基础，开展后续各项研究工作。图1展示了水系沉积物地球化学勘查采样点自动设计工作流程图。

图1 采样点位自动设计实现流程图Fig.1 Working flow chart of automatic design of sampling location

2.1.1 建立DEM

数字高程模型(Digital Elevation Model, DEM)是水系流域提取的基础，它通过将高程数据栅格化，可以很方便地计算栅格之间高程数据的关系，从而提取DEM中包含的地形、地貌、水文信息。

笔者采用ArcGIS 10.2软件，首先利用工作区内的高程点和等高线shape线文件数据创建TIN文件，再利用ArcToolbox工具箱内的Conversion Tools转换工具，将TIN数据转换为栅格，本文所采用的栅格为10 m×10 m。

2.1.2 加入已知水系

传统的水系自动提取方法中，直接通过DEM进行提取，发现自动提取的水系与实际水系的分布有较大出入，前人的研究也发现这一问题，并作出了相应的研究[10-12]，取得了一定的成果。

为了避免上述问题发生，我们提出了利用已知水系限制DEM的方法，很好地解决了这一问题。利用建立数字高程模型时的参数，对已知水系进行栅格化，然后通过栅格计算器，将没有水系的栅格设置为“0”，有水系的栅格设置为-100，然后将水系栅格与地形DEM进行栅格相加，从而得到已知水系限制的DEM。

2.1.3 水系自动提取

1) 水流方向计算。对比中心栅格与相邻栅格之间的高程，根据水往低处流的常理，假定水流向最大坡降的栅格流动，利用D8算法(D8算法的核心思想是：假定雨水降落在地形中某一个格子上，该格子的水流将会流向其周围8个格子中地形最低的格子中)，计算出DEM中所有栅格的流向。ArcGIS 10.2软件中，应用水文分析库下的流向确定命令,生成8方向水流流向，流向原理可参考文献[11]。

2)汇检测和洼地填补。流向确定后，所有的水流都向低处流动，然而工作区内的局部低洼点，周边栅格高程均高于局部低洼点，仅根据流向计算，所有的水系将终止于此栅格点。寻找低洼点的过程称为汇检测。汇的出现与实际现象相违背，因为低洼点在积满水后会向缺口的地方流动。因此找出洼地，然后填充洼地可以解决这一问题。

3) 栅格流量。填充洼地后的DEM，再计算各栅格的流向。假设每一个栅格因为大气降水产生流量为“1”的水量，按栅格流向向下游流动，每向下流动一个栅格，水量将增加“1”，而每一个栅格的流量便为所有流向此栅格的栅格个数。据此计算出每一个栅格点累积流量。

4) 设置水系流量阈值。实际上只有部分栅格会形成水系，而有的仅形成季节性冲沟。判定是否形成水系或者冲沟，最关键的问题是设置阈值。根据区内已知的水系资料，进行多次的试验可以给出与实际情况相符的阈值。为了区分不同流量的水系，还可以通过分类设置阈值来达到水系分类的目的。

5)水系提取。以设定的阈值对累积流量进行提取,大于阈值的栅格被定义为水系。

表1 水系特征点的分类

2.2 水系分类以及特征点提取

根据《地球化学普查规范》的要求[13]，每个采样点最好能控制0.25 km2面积的汇水域，由于本方法栅格化时采用的栅格为10 m×10 m，因此0.25 km2面积的汇水流域为2 500个栅格产生的流量。

为了使水系的出口有意义，采用了水系分类提取的办法。流量大于2 500为第一类；流量1 250至2 500为第二类；流量625至1 250为第三类；流量312至625为第四类(表1)。

显然此时第二类，流量为1 250至2 500，水系出口流量正好为2 500，此点为最符合《地球化学普查规范》的要求的水系沉积物采样点位。但是由于流量为2 500的点分布不均匀，对于山脊以及已知的大江大河内，此类点很难存在。因此为了符合《地球化学普查规范》要求中点位分布均匀的特点，必须在其他流量非2 500的点位上设计点位。

笔者提出的解决方法是：将各水系的出口点和长度大于500 m的第一类水系的中点提取出来，然后将这些提取出来的点按其所处的水系类别按表的规则为这些点赋予顺序号，以作为水系深积物设计点位筛选的关键字。

2.3 水系特征点内筛选设计点位

根据《地球化学普查规范》的要求[13]，累积流量为2 500的采样点位，控制的汇水面积正好是0.25 km2。但是此类点较少，仅筛选此类点很难满足采样密度要求。因此必须选取其他特征点作为设计采样点。

在同一个采样小格内，如果存在不同类型的特征点，则采取如下的处理方式进行筛选：①以顺序号为第一关键字升序排序，累积流量作为第二关键字降序排序的方式进行排序；②选取同一小格内的前两个特征点作为最终的采样设计点位。

2.4 设计点位自动编号

在1∶50 000水系沉积物地球化学勘查中，样品编号分为三步：

1)将图幅范围内不完整的方里格网坐标向图幅外延伸，取最邻近的整公里方里网为界，然后以1 km为间距纵向和横向上划分网格，并按从左到右，从上到下的顺序按数字编号。

2)将①划分好的网格以纵横方向500 m的间距将网格划分为四个小网格，左上角小网格编号为A，右上角小网格编号为B，左下角小网格编号为C，右下角小网格编号为D。

3)500 m×500 m的网格内样品再按从左到右，从上到下的顺序按从“1”开始的自然数编号，确定小格内样品顺序编号。

采样设计点最终编号为大格编号+小格编号+小格内样品顺序编号。自动编号的计算相对简单，只需要简单的数据换算和排序便可以完成。

假定图幅的四个边界坐标分别为left、right、top、bottom，待编号的样品设计点坐标为x、y，那么该点大格编号和小格编号通过运算可以得到。首先，通过excel函数求取四个边界的整公里坐标。

大格编号：(INT(top / 1000) - INT(y / 1000)) * (INT(right / 1000) - INT( left / 1000)) + INT(x/1000) - INT(left/1000)

小格编号：首先计算10 * INT( MOD( x , 1000) / 500)+ INT( MOD( y ,1000)/500)，根据计算结果选定小格编号，如果结果为“0”，则小格编号为C；如果结果为“1”，则小格编号为A；如果结果为“10”，则小格编号为D；如果结果为“11”，则小格编号为B。

小格内样品顺序编号主要是根据小格内样品从左到右，从上到下的顺序进行编号。首先，将计算得到的大格编号和小格编号按字符串连接起来，然后将其编号与其对应的样品坐标x、y组成一行，然后按excel的排序功能，将所有的设计点位进行排序，以大格编号和小格编号连接起来的字符串列为主要关键字，x为第一次要关键字，y为第二次要关键字。排序好之后的数据，通过判断当前行所对应的大格编号和小格编号连接起来的字符串与上一行的字符串是否相同，如果不同，则小格内样品顺序编号设置为“1”，如果相同，小格内样品顺序编号则为上一行小格内样品顺序编号加上“1”。

最终将大格编号，小格编号和小格内样品顺序编号按字符串连接起来便得到了样品编号。

3 试验区的应用及效果分析

3.1 试验区概况

此次研究选取了贵州省黔西南1∶50 000马岭幅作为试验区。

马岭幅位处南盘江-右江前陆盆地大地构造(相)之兴义泥凼-贞丰者相台缘-斜坡亚相。由师宗-弥勒断裂带和紫云-六盘水断裂带在贵州省内围限区域，属滇黔桂“金三角”的组成部分，是由扬子被动边缘碳酸盐台地演化而成的一个中晚三叠纪周缘前陆盆地。测区地层为上古生界至中生界，其中以三叠系中上统的陆源碎屑复理石最引人瞩目。测区构造变形组合样式具显著“条”、“块”镶嵌特色，即北东向，北西向强变形带(“条”—紧闭形褶皱及逆冲断裂)与弱变形块(宽缓短轴型褶皱区)镶嵌。

马岭幅内目前尚未发现金矿或金矿点，但是图幅外分布着许多金矿或者金矿点，如图幅北侧有泥堡金矿，楼下金矿，东侧有戈塘金矿，南侧有雄武金矿等。在区域地球化学金异常图中，马岭幅内不仅有幅值较小的金异常，而且区内北东向，北西向和近东西向断裂发育，区内出露二叠系龙潭组地层。因此，在马岭幅开展1∶50 000水系沉积物地球化学勘查，对寻找以陆源硅质碎屑岩为主要容矿岩石的微细浸染型金矿床有一定的潜力。

图2 自动提取水系结果对比图Fig.2 A contrast diagram of water system extraction

3.2 采样点位自动布置效果

3.2.1 不同方法的水系提取效果对比

笔者在对水系自动提取的处理上，相较于传统的自然水系提取，采用了已知水系限制后再提取水系的方案，取得的结果与实际水系重合度更高，更好地保证了后期采样点位置的精度。

没有已知水系限制下的常规传统方法提取的水系，在地势起伏较大的地段,与已知水系分布较为吻合；而在宽阔山谷中，河道位置偏差较大，尤其在已知河道转弯处，吻合度很低 ，提取的水系不能用作采样点位布置的参考。

应用本文提出的已知水系限制条件下自动提取的水系，不仅在地势起伏较大的地段吻合，而且在宽阔山谷中，水系弯道段都能有很好地吻合，大大提高了自动提取水系的精度，使得后期点位设计的位置更精准。

图3 采样点位自动布置结果图Fig.3 Distribution map of sampling point

3.2.2 采样点自动布置结果

图3为使用本文提出的方法自动布置的采样点位(因为一个图幅较大，为便于展示，仅提取了一个流域)，采样点位布置均匀，采样点均位于水系或者是季节性水系上，很好地控制了采样点上游的汇水域，采样点编号正确，根据水系分布密度不同，采样小格的采样点最多为2个，空白小格数较少，可以满足水系沉积物测量采样点位设计的基本要求。但是，图3中尚有不少设计点位的合理性不够，如3级水系中设计点位太多等，需要今后开展进一步的研究工作，完善本方法技术。

4 结论

笔者结合GIS和水文科学相关知识，为水系沉积物地球化学勘查采样点位布置提供了一种自动、快捷的新方法。以1∶50 000电子版地形图为基础数据，利用ArcGIS 10.2软件平台，通过已知水系的限制，自动提取出与实际情况相吻合的水系。根据流量将水系进行分类，并提取水系的出水口和流量大于2 500 m且长度大于500 m水系的中点。筛选出符合《地球化学普查规范》中采样点要求的特征点作为最终采样点。通过采样点坐标换算，自动进行编号，实现了水系沉积物勘查采样点位自动布置。

新方法是多学科多领域交叉融合的结果，极大地提高了工作效率，大大减轻了从业人员劳动量。同时，克服了不同经验工作人员布置点位不同的缺点。

由于水系自动提取方法自身的缺陷，在喀斯特地区和平顶山上，尤其是岩溶漏斗和“断头河”较多的地区，该方法误差较大。同时，由于选择条件设置不够全面，筛选出来的设计采样点位还有不少不符合《地球化学普查规范》的要求，需要人工干预调整。

这里仅提出了一种思路，且是初步尝试，具有了一个很好的开端，以后还需要根据野外实际情况，提取河道更多的特征点。同时，按照《地球化学普查规范》中的布点原则，设置更多的选择条件，最终达到筛选出完全符合实际要求的水系沉积物测量设计采样点。