CN105005649A - 一种进行潜在地下水分布制图的方法 - Google Patents

Abstract

本发明公开了一种进行潜在地下水分布制图的方法，首先利用DEM数据，提取各种地下水制图的特征因子；根据各因子之间相互影响关系，建立累积影响因子矩阵，用于决定各特征因子之间的权值分配；将图形统一转换为256级灰度图，根据特征因子和权值分配，采用多准则分析方法，利用ArcGIS中的加权叠加工具，进行地下水潜在值计算并制作初步结果图；对初步结果图进行均值滤波，以减少、过滤噪声信息，然后根据自然断点法进行分级分类，将初步结果图分为四个等级。本发明的有益效果是采用比较容易获取的DEM数据作为单一数据源，提取潜在地下水分布相关特征因子，进行潜在地下水分布建模。验证结果显示模型有比较高的成功率。

Description

一种进行潜在地下水分布制图的方法

技术领域

本发明属于水资源勘探技术领域，涉及一种进行潜在地下水分布制图的方法。

背景技术

地下水是最为重要的自然资源之一，可为人们提供经济、可靠的饮用水。目前，地下水供水大概占年供水量34％(N.S.Magesh 2012)，是一种非常重要的淡水资源。因此对于地下水资源进行探索评估是地下水可持续管理发展的需求。特别是在干旱地区，进行地下水的准确探测以进行开发利用，显得尤为必要。

近年来，众多学者发现利用“多准则决策技术”(Multi-criteria decisionmaking technique)在进行潜在地下水分布评估中是一种非常有效的工具。如Giuseppe Mendicino等研究后提出地中海气候条件下的地下水资源指数(Giuseppe Mendicino 2008)，Ratnakar Dhakate等则结合地貌学和地球物理方法来定位潜在地下水区域(Ratnakar Dhakate 2008)，Manika Gupta等结合GIS、RS和相关野外技术，联合进行潜在地下水的定位实验(Manika Gupta 2010)，Hyun-Joo Oh等利用GIS和概率模型技术评估潜在地下水含量的研究(Hyun-JooOh 2011)，Saro Lee等利用GIS和证据权方法(weights-of-evidence method)进行了区域潜在地下水产出的制图研究(Lee,Kim et al.2012)。Samy IsmailElmahdy和Mohamed Mostafa Mohamed结合GIS和RS方法，利用DEM和LandsatTM影像进行了潜在地下水区域的划分(Mohamed 2014)，D.Davoodi Moghaddam等利用二元统计模型结合GIS进行了潜在地下水的制图研究(Moghaddam,Rezaeiet al.2015)，等。

尽管多准则决策技术被广泛应用于潜在地下水分布描述中，但是该技术需要尽可能多的数据和信息，以构成多种准则以综合进行评判。这在数据稀少地区和干旱地区的研究则会面临问题。由于缺少相关数据，多准则决策技术需要的多个参数条件难以满足，使得研究难以为继。

DEM(Digital Elevation Model，数字高程模型)是地表形态的数字化表达。由于GIS技术的发展以及越来越多的高质量的DEM数据的获取，使得DEM的应用范围越来越广，如水力水文、地形地貌、生态环境等众多领域的研究中(Thompson,Bell et al.2001)。针对DEM的研究和应用众多，现有众多文献提出了如何从DEM数据中提取众多信息，如涉及地下水研究中的地形因素、水文因素、地质因素等，都可以从DEM数据中提取得到，这使得利用DEM进行信息提取以进行潜在地下水分布制图成为可能。

发明内容

本发明的目的在于提供一种进行潜在地下水制图的方法，解决了现有的潜在地下水制图需要尽可能多的数据和信息，以构成多种准则以综合进行评判，而这在数据稀少地区和干旱地区由于缺少相关数据，多准则决策技术需要的多个参数条件难以满足，使得研究难以为继的问题。

本发明所采用的技术方案是按照以下步骤进行：

步骤1：首先利用DEM数据，提取各种地下水制图的特征因子；

步骤2：根据各因子之间相互影响关系，建立累积影响因子矩阵，用于决定各特征因子之间的权值分配；

步骤3：潜在地下水制图：

将图形统一转换为256级灰度图，根据特征因子和权值分配，采用多准则分析方法，利用ArcGIS中的加权叠加工具，按照下式进行地下水潜在值计算并制作初步结果图：

GPV＝D_l×31％+TWI×23％+Relief×8％+D_D×15％+CI×23％；

步骤4：对初步结果图进行均值滤波，以减少、过滤噪声信息，然后根据自然断点法进行分级分类，将初步结果图分为四个等级：非常好潜在区域、好的潜在区域、中等潜在区域、差的潜在区域。

本发明的有益效果是采用比较容易获取的DEM数据作为单一数据源，提取潜在地下水分布相关特征因子，进行潜在地下水分布建模。验证结果显示模型有比较高的成功率。

附图说明

图1是研究区线性地质断裂特征提取结果及其风向玫瑰图；

图2是线性地质断裂分布密度图；

图3是研究区坡度分布图；

图4是研究区坡度累计百分比分布曲线；

图5是CI计算示意图；

图6是CI＝-90°示意图；

图7是CI＝+90°示意图；

图8是CI＝0°示意图；

图9是研究区不同区域矩形大小VS地形起伏散点图；

图10是潜在地下水预测值的灰度直方图及自然断点法的断点示意图；

图11是研究区潜在地下水分布图；

图12是潜在地下水成功比例曲线示意图。

具体实施方式

下面结合具体实施方式对本发明进行详细说明。

本发明利用DEM数据提取与潜在地下水分布密切相关特征因子信息，然后利用“多准则决策技术”(Multi-criteria decision making technique)，给各因素进行合理权重分配后进行组合叠加得到最终结果。本发明的具体步骤如下：

步骤1：首先利用DEM数据，提取各种地下水制图的特征因子；

本发明中，提取5种因子，即地质断裂特征、地形湿度指数、地形收敛指数、河网特征、地形起伏特征，以进行以进行潜在地下水分布建模。

地质断裂特征

地表的线性地质断裂特征是地底结构特征(如硬岩地区的断裂、裂缝)的反映，是地下水径流的通道，因而是重要的地下水存在指示标志。一般说来，断裂构造密度越大，蕴含地下水的可能性越大。线性地质断裂在地下水的出现和运动中具有非常重要的角色，因此被选取用来作为潜在地下水分布的特征因子之一。

从DEM中提取线性地质断裂特征，Guth(1999,2001)提出一种算法用以从DEM中提取地形断裂特征。该算法需要四个输入参数，即点间隔、区域大小、平坦去除阈值和长度倍数。点间隔是利用DEM数据进行纹理计算时的点之间的距离(单位：米)，点间隔越大，提取出的线性地质断裂特征越少；区域大小(单位：米)是利用DEM数据提取信息时进行运算的区域块的尺寸，区域尺寸过小的话会得到不合理的结果；平坦去除阈值是指大于该阈值的点将不会作为地形断裂特征进行绘制，所以该值越大，绘制的点将越多，但这些点也可能属于随机噪声点。所以合适阈值的选取很关键，一般要根据DEM数据的比例以及研究地区的地形起伏情况进行选择；第四个参数是长度倍数，是提取出的线性纹理矢量的比例因子，该值越大，提取出的线性地质断裂将会越长。根据该算法，进行研究区的线性地质断裂特征的提取结果如图1所示。

根据图1可以看出，该地区的线性地质断裂特征在西部和中部相对集中，主要成西北—东南和正北—正南走向。提取出研究区线性地质断裂后，可以进一步进行其分布密度制图。线性特征的密度可以用线性特征的总长度除以该地区总面积得到，如式(1)所示：

$D_1={\Σ}_{i=1}^n{L}_i/A---\left(1\right)$

其中L_i代表第i条线性特征的长度，A表示研究区的总面积，D_l即为线性特征密度。据此可以计算得到其密度分布图如图2所示：

地形湿度指数

由于地形在水文过程中的重要作用，地形湿度指数(topographic wetnessindex，TWI)的概念被提出，用于定量描述地形因素对于水文过程的控制。由于能模拟流域内土壤水分的干湿状况，在流域的土壤及分布式水文模型等研究中具有重要的意义。TWI能反映潜在地下水由于地形因素而导致的渗漏，故被选取用来进行潜在地下水分布建模。

TWI是上坡汇水面积和坡度的函数，其定义如下式(2)所示：

TWI＝ln(α/tanβ)     (2)

式中，α为流经地表某点的单位等高线长度上的特定汇水面积(specificcatchment area,SCA,m².m^-1)，β为该点的坡度角。对于栅格形式的DEM数据，α为网格单元汇水面积与DEM栅格尺寸的比值。TWI越大，反映出该点坡度越平缓，上坡汇水面积也越大。在地形起伏较为明显的山区，该式能够得到比较理想的结果。但是在地形平坦地区，该式则差强人意，有时会得到明显不合理的结果。因此，众多学者对其进行了改进，其中，Bohner等提出的改进算法(Bohner&Selige,2006)效果显著。该算法针对汇水面积进行了改进，与传统TWI算法单纯计算α不同，该算法对小于一定值的特定特定汇水面积(SCA)的计算结果α进行了迭代运算直至其大于给定值，其定义如下式(3)所示，用SCA_M代替α，得到最后的地形湿度指数计算公式如式(4)所示。

$\begin{array}{ccc}{\mathrm{SCA}}_M=SCAmax{\left(\frac{1}{15}\right)}^{\mathrm{\&beta;}\exp \left({15}^{\mathrm{\&beta;}}\right)}& For& SCA<SCAmax{\left(\frac{1}{15}\right)}^{\mathrm{\&beta;}\exp \left({15}^{\mathrm{\&beta;}}\right)}\end{array}---\left(3\right)$

$TWI=ln\left(\frac{{\mathrm{SCA}}_M}{tan\mathrm{\&beta;}}\right)---\left(4\right)$

如图3所示为研究区的坡度分布图，可以看出坡度较大的地区分布面积较小，主要集中在西部高地，其它大部分地区坡度都较为平坦；图4为研究区坡度累计百分比分布曲线，可以看出有95％的坡度都小于6％。可见研究地区属于平坦地区，因此采用式(4)进行地形湿度指数的计算。

地形收敛指数

地形收敛指数(Convergence Index，CI)，是用来描述地表径流的汇聚(收敛)区域和发散区域)(Kiss 2004(N.Thommeret,2009#3)),因此可以用来进行地下水分布建模。CI可以基于从DEM中提取的坡向数据进行计算。通过计算相邻像元坡向和指向中心像元的夹角的平均值，减去90°来计算得到。

如图5所示为CI计算示意图，中心点(黑色圆点)的收敛指数计算公式如下式(5)计算得到：

其中θ_i为第i个中心像元点的相邻像元坡向和指向中心像元的夹角，CI即为收敛指数。收敛指数的取值范围为[-90°，+90°]，正值代表径流发散区域，负值代表径流收敛区域，0则代表平地。如图6所示，θ为0°，CI为-90°的情况，图7是θ为180°，CI为90°的情况，图8是θ为90°，CI为0°的情况。如此，根据如上计算方法，首先从DEM数据中提取坡向数据，再根据式(5)计算得到其收敛指数。

河网特征

从DEM数据中提取河网特征的研究很多，算法也很成熟(Kiss 2004)。在众多算法中，D8算法应用广泛，并且被ArcGIS软件采用(Mohamed 2014)。本发明中采用D8算法，运用ArcGIS软件进行计算提取得到区域的河网分布图。

根据提取出的河网分布，可以进一步计算研究区的河网分布密度。河网密度可以反映该地区的渗透度，因而被选取用来进行潜在地下水分布制图。低密度的河网地区表明该地区地表可渗透性强，有利于地下水的分布；反之高密度的河网地区则表明该地区地表可渗透行低，不利于地下水的分布。本研究中，利用ArcGIS中的线密度分析工具进行河网密度分布制图。河网密度可以用D_D(Drainage Density)来表示。

地形起伏特征

地形起伏(Relief)在决定降水的渗透速率，径流汇聚、传输和损耗区域中有着重要作用(Mohamed 2014)，因此可用于潜在地下水分布制图中。地形起伏特征可以用指定大小的区域内(通常为一矩形区域)的最大高程差来反映(如式6所示)：

Relief＝MaxElevation-MinElevation  (6)

MaxElevation为指定区域内最大高程值，MinElevation则为该指定区域内的最小高程值，Reilief即为指定区域内的地形起伏度。合适的区域矩形大小的选择很重要。在本研究中，如图9所示的地形起伏(Relief)vs区域矩形大小(Region Box)的散点分布图，可见随着区域矩形尺寸的增加，地形起伏也在变大。在区域矩形大小为600米之前，地形起伏随着区域尺寸的增加而快速变大；在区域矩形大小为600之后，地形起伏随着区域矩形的增加而变大的趋势则较为平缓。因此本研究选取区域矩形大小为600米进行研究区的地形起伏制图。

步骤2：根据各因子之间相互影响关系，建立累积影响因子矩阵，并据此计算各因子权重分配；根据权值分配进行叠加分析，得到初步结果图；

提取了潜在地下水分布的特征信息之后，需要进行各特征因子的权值分配。传统的权值分配，由于数据充足，可以将钻井数据分为2部分，一部分作为训练样本数据，一部分作为验证数据，来获取最佳权值分配。同样，该方法在数据比较稀缺地区无法得到很好的应用。本研究中，根据各特征因子之间的相互影响关系，建立对应的“累积影响矩阵”，用于决定各特征因子之间的权值分配。

如表1所示，为步骤一所提到的5个特征因子所组成5×5的矩阵，以及权值累积结果(Σ列)和权值分配结果(权值列)。权值分配按照如下规则进行：对于5×5的矩阵中的每个特征因子，如果某特征因子对另外的特征因子有影响，则在对应的交叉矩阵格中赋值为1，反之为0。对角线上的元素统一为1(即特征因子对自身有影响)。以“地质线性断裂特征”为例，由于“地质线性断裂特征”会影响到相应的沟谷、山脊分布，因此对于矩阵中的特征因子“地形湿度指数(TWI)”、“河网分布”和“地形收敛指数(CI)”都有影响，但对“地形起伏(Relief)”无影响，故在表1中的各格中赋值依次为1、1、0、1、1，累计值(Σ)为4。如此分别计算出各特征因子的累积影响因子，各特征因子的权值分配可由该特征因子的累积影响因子除以所有特征因子的累积影响因子得到。

表1：累积影响矩阵及权值分配

步骤3：潜在地下水制图；将图形统一转换为256级灰度图；根据各特征因子和权值分配，采用多准则分析方法，利用ArcGIS中的加权叠加工具，按照下式进行地下水潜在值(Groundwater Pontential Value，GPV)计算，并据此值得到结果图。

GPV＝D_l×31％+TWI×23％+Relief×8％+D_D×15％+CI×23％。

对结果图进行均值滤波，以去除噪声点；并根据“自然断点”法进行分类，将整个区域分为4个等级：非常好潜在区域、好的潜在区域、中等潜在区域、差的潜在区域。同时注意根据DEM值以及TWI计算结果提取出非潜在区域(如图10所示，在本研究区中主要为水体区域)。用来进行分类的直方图及最终结果图如图11所示。根据图11可知，研究区潜在地下水分布中，非潜在区域、非常好潜在区域、好的潜在区域、中等潜在区域、差的潜在区域分别占研究区总面积的2.6％，19.8％，39.4％，29％和9.2％。其中非潜在区域主要是水体区域，非常好和好的潜在区域基本分布在研究区西部，也即地质线性断裂密度高区域；而差的潜在区域主要分布在研究区东部。

结果验证

根据该地区的钻井数据，将钻井位置和潜在地下水分布地图相叠加，从潜在地下水地图的最高值开始计算地下水出现的累积百分比。以累积潜在地下水地图灰度值为x轴，以累积地下水出现百分比为y轴，将计算得到的结果进行绘制得到成功比例曲线如图12所示。

根据接受者操作特性曲线(receiver operating characteristic，ROC)理论，计算曲线以下面积(area under curve，AUC)占整个坐标矩形范围的比例，来进行结果验证。经计算，AUC为0.79，即本研究所建立模型计算得到的潜在地下水地图成功预测的准确率为79％，预测结果可以接受。

本发明的优点还在于：本研究针对地下水研究中数据缺少的状况，采用比较容易获取的DEM数据作为单一数据源，提取潜在地下水分布相关特征因子，进行潜在地下水分布建模。验证结果显示模型有比较高的成功率。该研究可为数据稀少地区的地下水研究开辟一条新的思路，同时也拓展了DEM数据的应用范围。

以上所述仅是对本发明的较佳实施方式而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施方式所做的任何简单修改，等同变化与修饰，均属于本发明技术方案的范围内。

Claims (1)

1.一种进行潜在地下水分布制图的方法，其特征在于按照以下步骤进行：

步骤1：首先利用DEM数据，提取各种地下水制图的特征因子；

步骤2：根据各因子之间相互影响关系，建立累积影响因子矩阵，用于决定各特征因子之间的权值分配；

步骤3：潜在地下水制图；

将图形统一转换为256级灰度图，根据特征因子和权值分配，采用多准则分析方法，利用ArcGIS中的加权叠加工具，按照下式进行地下水潜在值计算并制作初步结果图：

GPV＝D_l×31％+TWI×23％+Relief×8％+D_D×15％+CI×23％；

步骤4：对初步结果图进行均值滤波，以减少、过滤噪声信息，然后根据自然断点法进行分级分类，将初步结果图分为四个等级：非常好潜在区域、好的潜在区域、中等潜在区域、差的潜在区域。