CN107479042B - 一种表层岩溶带空间蓄水能力的估算方法 - Google Patents

Abstract

本发明公开了一种表层岩溶带空间蓄水能力的估算方法。括下述步骤：a、利用探地雷达获得表层岩溶带雷达信息，同时记录采样点位的空间信息；b、对表层岩溶带雷达信息进行滤波处理，利用层位追踪的方法提取出表层岩溶带雷达信息中的深度信息；c、根据空间信息和深度信息两者存在的空间体积差值获取表层岩溶带的体积；d、获取表层岩溶带的裂隙率，计算出表层岩溶带的裂隙体积；核算表层岩溶带裂隙中的被土壤充填体积与未被土壤充填体积所占比例，计算被土壤充填裂隙的持水体积，被土壤充填裂隙的持水体积与未被土壤充填体积之和即为表层岩溶带的蓄水体积，即蓄水能力。本发明能够有效估算出表层岩溶带空间蓄水能力，计算方法更加简单。

Description

一种表层岩溶带空间蓄水能力的估算方法

技术领域

本发明涉及物探领域应用于空间体积估算，特别是一种表层岩溶带空间蓄水能力的估算方法。

背景技术

表层岩溶带是喀斯特地区重要的蓄水空间，起到储存、滞留、调蓄水分的关键作用，估算其蓄水能力对于认知喀斯特区地下结构、水循环、物质循环研究具有重要的理论意义，通过该方法定量评估、开发与表层岩溶带的水资源，解决喀斯特地区工程性缺水的现状也具有重要的现实意义。目前对于表层岩溶带蓄水能力的定量评估国内外仍处于研究空白。

发明内容

本发明的目的在于，提供一种表层岩溶带空间蓄水能力的估算方法。本发明能够有效估算出表层岩溶带空间蓄水能力，计算方法更加简单。

本发明的技术方案：一种表层岩溶带空间蓄水能力的估算方法，包括下述步骤：

a、利用探地雷达获得表层岩溶带雷达信息，同时记录采样点位的空间信息；

b、对所述的表层岩溶带雷达信息进行滤波处理，利用层位追踪的方法提取出表层岩溶带雷达信息中的深度信息；

c、根据空间信息和深度信息两者存在的空间体积差值获取表层岩溶带的体积；

d、获取表层岩溶带的裂隙率，计算出表层岩溶带的裂隙体积；核算表层岩溶带裂隙中的被土壤充填体积与未被土壤充填体积所占比例，计算被土壤充填裂隙的持水体积，所述的被土壤充填裂隙的持水体积与未被土壤充填体积之和即为表层岩溶带的蓄水体积；

e、步骤d所述的表层岩溶带的蓄水体积即为表层岩溶带空间蓄水能力。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤a中，所述的利用探地雷达获得表层岩溶带雷达信息，包括下述步骤：

a1、设置数据采集窗口；

a2、通过探地雷达的发射天线与接收天线探测，在数据采集窗口反映出原始雷达影像；

a3、对比原始雷达影像所反映出的表层岩溶带的剖面信息，确定表层岩溶带的介电常数及层位划分原则，得到表层岩溶带深度信息。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤b中，所述的进行滤波处理，是利用Reflex处理软件进行滤波处理，具体包括下述步骤：

b1、导入数据，设置表层岩溶带的介电常数；

b2、一维滤波/去直流漂移，除去零点漂移，设置两个参数，所述的两个参数分别为1time.ns和2time.ns；所述的1time.ns数值为时间窗口2/3处的值；所述的2time.ns数值为从最后样点倒数第5个样点处的数值；

滤波公式为：

其中，ω₁(t)是处理前第1道扫描回波与ω₁′(t)是处理后第1道扫描回波，M是扫描回波的总道数，ω_l(t)是处理前第l道扫描回波；

b3、静校正/移动开始时间，通过分析得到的第一个波峰处为直达波时间，把直达的开始位置移到零点处；

b4、增益/能量衰减，对探地雷达信号放大，选择所需增益放大的比例因子；

b5、二维滤波/抽取平均值

设置输入开始时间start time和结束时间endtime；滤波公式为：

计算，其中Wm,n(t)为多次测量接收到的信号，m为各点测量序号，n为重复次数，M为扫描回波的道数；

b6、一维滤波/巴特沃斯带通滤波，设置低切频率与高切频率，得到有效波的图像；滤波公式为：

其中，n为滤波器阶数，ω_c为截止频率，ω_p为通频带边缘频率；H(ω)为振幅；

b7、二维滤波/滑动平均；滤波公式为：

其中n为滤波器阶数，ω_c为截止频率，ω_p为通频带边缘频率；H(ω)为振幅。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤c中，所述的空间体积差值，是步骤a所述的空间信息中的地面高程信息与步骤b所述的深度信息存在的空间体积差值。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤c中，所述的空间体积差值，按下述步骤计算：

c1、选取基准面作为计算表层岩溶带的体积的底界面，并在探地雷达获得表层岩溶带雷达信息后，运用空间插值，获得表层岩溶带的表层岩溶带下界；

c2、分别计算基准面与地面、基准面与表层岩溶带下界所构成的空间结构的体积，两者之差为空间体积差值，即表层岩溶带的体积。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤c1中，所述的选取基准面，是选取低于表层岩溶带下界的基准面。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤d中，所述的被土壤充填裂隙的持水体积为土壤田间持水度与土壤本身含水量的差值。

前述的表层岩溶带空间蓄水能力的估算方法所述的步骤a中，所述的记录采样点位的空间信息，是使用GPS记录采样点位的空间信息。

有益效果

与现有技术相比，本发明利用空间体积差值计算表层岩溶带的体积，并结合表层岩溶带裂隙的蓄水能力，计算得到表层岩溶带的蓄水能力，其方法简单。本发明将表层岩溶带裂隙分为未被土壤充填裂隙与被土壤充填裂隙，在估算表层岩溶带裂隙的蓄水能力时，将未被土壤充填裂隙按最大蓄满原则得到未被土壤充填裂隙的蓄水能力，再计算被土壤充填裂隙的持水体积，两者(即被土壤充填裂隙的蓄水能力、未被土壤充填裂隙的持水体积)之和即估算得表层岩溶带裂隙的蓄水能力，即为表层岩溶带的蓄水能力；通过该方式，计算得到的表层岩溶带的蓄水能力更加精确。

本发明采用探地雷达对表层岩溶带的深度探测替代了传统的典型剖面调查、钻孔调查，因而本发明具有速度快、分辨率高、操作方便、探测成本低、无破坏性和探测成本低的优点；不仅如此，本发明采用探地雷达探测，使得本发明能对表层岩溶带的深度信息进行定量化分析，并能够结合RTK空间定位获得高精度采样点位的点面信息，进而提高探测精度。

本发明通过探地雷达探测，使得本发明对于表层岩溶带的识别适用于任何地区，本发明利用RTK点位信息(即空间信息)与表层岩溶带的深度信息结合，运用空插值分析，充分认识了表层岩溶带在空间分布的特征，继而能够对整个区域内进行分析，而不再局限于某一剖面，很好地估算出了表层岩溶带的空间体积信息，进而计算其蓄水能力更方便。本发明通过表层岩溶带典型剖面调查与计算，以得到表层岩溶带的裂隙的各项信息(包括表层岩溶带的裂隙率、裂隙体积、被土壤充填裂隙的持水体积、未被土壤充填体积等)，克服了利用探地雷达探测带来的因表层岩溶带的异质性强、介电常数的确定难度大，对于表层岩溶带层位的提取及解释较为困难的问题，进一步使得表层岩溶带的蓄水能力的计算更加方便和精确。

为了证明本发明有益效果，申请进行如下试验：

通过50MHz的探地雷达对野外样地进行探测，获得原始图像及GPS空间点位信息。

采集的数据结果通过Mala公司提供的REFLEXW软件对原始图像进行滤波处理，主要通过直流漂移、静校正/移动初始时间、增益/能量衰减、二维滤波/滑动平均、1维滤波/带通滤波、抽取平均值等步骤，获得表层岩溶带滤波结果图像，然后通过相位自动追踪识别，提取出表层岩溶带的厚度结果(流程如图3所示)。本文选择剖面11来阐述雷达解译流程。

表层岩溶带的介电常数为0.1m/ns。

具体流程如下：

1.去直流漂移(一维滤波)

选择时间一为总时窗的三分二处，时间二选择为离雷达最后样点5个样点处，因为此地区总样点都约在612左右，为了统一所有剖面的一致性，所以全部剖面两个时间值分别选择为612、408，经过对波形的失真处理修正后，发现雷达波图像的信噪比增强，更能清楚的识别深层信息，得到结果图像如图4所示。

2.静校正(移动开始时间)

静校正处理后的雷达图像如图5所示。

3.增益

因为此次探测的深度主要是在30m以内，为了更加清楚的识别表层岩溶带的信号，并且在石灰岩地区所以Scaling value选择1左右(越大放大的倍数就越大)，得到增益处理后的图像如图6所示。

4.二维滤波(抽取平均值)

该滤波作用主要用于选定的道数，对其结果图像中每个时间段中选定的道数进行抽取平均，其执行方式为滑动背景的去除。考虑此次探测过程中剖面的的长度不同，所以选定的道数为(参数平均道average traces)为100，则根据计算原则选择左右两测点各50个点考虑；同时还限制了滤波作用的时间范围，并且选定参数开始时间与结束时间，对于数据结果对应道数的其它部分不进行滤波处理，作为缺省值每道的开始时间设为数据的起点，结速时间设为每道最大时间处，可以抑制水平一致能量，其效果主要表现在存在差异的信号部分。其处理结果如图7所示。

5.巴特沃斯带通滤波

在探测表层岩溶带厚度的过程中，选择的天线为50MHz，所以低切频率选择18MHz，高切频率为60MHz，其余的频谱都设定为0，通过处理得到的图像如图8所示。

6.二维滤波(滑动平均)

其参数平均道数(average traces)一般取值为3，同时还可以限制滤波作用的时间范围，即参数开始时间与结束时间，因为喀斯特的高异值性，所以时间范围一般选取缺省值，对于超过此部分的数据不进行滤波作用，得到的雷达图像如图9所示。

通过以上处理得到清晰的雷达图像结果，采用自动相位识别进行层位追踪，识别出表层岩溶带层厚度。得到11号剖面表层岩溶带厚度层位图(如图10所示)，结果显示，此次总计剖面35条，通过空间插值获得表层岩溶带的深度面空间图(如图11所示)。

结论：通过基准面的选择，计算得到该区域内表层岩溶带的体积为2028542.2立方米，由典型剖面裂隙率的统计分析发现，其总裂隙率为0.083，通过计算该区内总的表层岩溶带裂隙体积为168369.003立方米，计算出表层岩溶带的蓄水体积最大可以达到17670立方米，表明表层岩溶带的蓄水能力相对较强，在此研究区内能产生剖面径流需要最小降雨量为52毫米，与中国科学院地球化学研究于普定喀斯特生态系统观测研究站监测结果为60毫米相差不大，主要是因为此过程地下流失等其它影响因素不能确定，得到结果小于实际监测结果，基本与实际情况相似。

附图说明

图1是本发明的估算流程图；

图2是表层岩溶带的结构示意图；

图3是表层岩溶带波滤波流程图；

图4是去直流流漂移处理后的雷达图像；

图5是静校正处理后的雷达图像；

图6是增益处理后的雷达图像；

图7是二维滤波处理后的雷达图像；

图8是巴特沃斯带通滤波处理后的雷达图像；

图9是滑动平均处理后的雷达图像；

图10是11号剖面表层岩溶带厚度图；

图11是表层岩溶带的空间发育图。

具体实施方式

下面结合附图和实施例对本发明作进一步的说明，但并不作为对本发明限制的依据。

实施例1。一种表层岩溶带空间蓄水能力的估算方法，如图1所示，包括下述步骤：

a、利用探地雷达对表层岩溶带的深度探测得到表层岩溶带雷达信息，同时记录采样点位的空间信息；

b、对所述的表层岩溶带雷达信息进行滤波处理，利用层位追踪的方法提取出表层岩溶带雷达信息中的深度信息；

c、根据步骤a所述的空间信息和步骤b所述的深度信息，利用空间体积差值获取表层岩溶带的体积；

d、获取岩石的裂隙率，计算出表层岩溶带的裂隙率及裂隙体积；核算表层岩溶带裂隙中的被土壤充填体积与未被土壤充填体积(根据最大蓄满原则，未被土壤充填裂隙全部蓄水)所占比例，及核算被土壤充填裂隙的持水体积，所述的被土壤充填裂隙的持水体积与未被土壤充填体积之和即为表层岩溶带的蓄水体积。通过野外实际调查对典型的表层岩溶带剖面进行分析获得岩石的裂隙率，计算出表层岩溶带的裂隙体积(根据裂隙率与表层岩溶带的体积间的积即为裂隙体积)；分析表层岩溶带裂隙中被土壤充填裂隙与未被土壤充填裂隙的体积所占的比例(分析的方法是，通过典型剖面调查计算获得)。所述的被土壤充填体积与未被土壤充填体积所占比例采用出露的典型剖面分析计算获得。

e、步骤d所述的表层岩溶带的蓄水体积即为表层岩溶带空间蓄水能力。

前述的步骤a中，所述的利用探地雷达对表层岩溶带的深度探测得到表层岩溶带雷达信息，包括下述步骤：

a1、设置数据采集窗口；数据采集窗口设置包括：探地雷达选择使用的天线频率、触发数据采集方式(距离、时间或点测)及直达波零点；

a2、通过探地雷达的发射天线与接收天线探测，在数据采集窗口反映出原始雷达影像；

a3、对比原始雷达影像所反映出的表层岩溶带的剖面信息，确定表层岩溶带的介电常数及层位划分原则，得到表层岩溶带深度信息。

前述的步骤b中，所述的进行滤波处理，是利用Reflex处理软件(雷达数据处理的一个专用软件)进行滤波处理，具体包括下述步骤：

b1、导入数据，设置表层岩溶带的介电常数；所述的数据，是采集到的表层岩溶带雷达信息；

b2、一维滤波/去直流漂移，除去零点漂移，设置两个参数，所述的两个参数分别为1time.ns和2time.ns；所述的1time.ns数值约为时间窗口2/3处的值；所述的2time.ns数值为从最后样点倒数第5个样点处的数值；上述的两个参数根据探测深度、雷达频率的选择而不同。

滤波公式为：

其中，ω₁(t)是处理前第1道扫描回波与ω₁′(t)是处理后第1道扫描回波，M是扫描回波的总道数，ω_l(t)是处理前第l道扫描回波；

b3、静校正/移动开始时间，通过分析得到的第一个波峰处为直达波时间，把直达的开始位置移到零点处；

b4、增益/能量衰减，对探地雷达信号放大，选择所需增益放大的比例因子；比例因子因深度的不同而不同；

b5、二维滤波/抽取平均值

设置输入开始时间start time和结束时间endtime；滤波公式为：

计算，其中Wm,n(t)为多次测量接收到的信号，m为各点测量序号，n为重复次数，M为扫描回波的道数；

b6、一维滤波/巴特沃斯带通滤波，设置低切频率与高切频率，得到有效波的图像；滤波公式为：

其中，n为滤波器阶数，ω_c为截止频率，ω_p为通频带边缘频率；H(ω)为振幅；高低切频率的选择会因为雷达天线频率的选择而不同,对应的雷达频率会有固定的高低切频率值；

b7、二维滤波/滑动平均；滤波公式为：

其中n为滤波器阶数，ω_c为截止频率，ω_p为通频带边缘频率；H(ω)为振幅。

前述的步骤c中，所述的空间体积差值，是步骤a所述的空间信息中的地面高程信息与步骤b所述的深度信息存在的空间体积差值。

前述的步骤c中，所述的空间体积差值，按下述步骤计算：

c1、选取基准面作为计算表层岩溶带的体积的底界面，并在探地雷达获得表层岩溶带雷达信息后，运用空间插值，获得表层岩溶带的表层岩溶带下界；

c2、分别计算基准面与地面、基准面与表层岩溶带下界所构成的空间结构的体积，两者之差为空间体积差值，即表层岩溶带的体积。

前述的步骤c1中，所述的选取基准面，是选取低于表层岩溶带下界的基准面。表层岩溶带的结构如图2所示。

前述的步骤d中，所述的被土壤充填裂隙的持水体积为土壤田间持水度与土壤本身含水量的差值。土壤田间持水度与土壤本身含水量通过实际测试得到。

前述的步骤a中，所述的记录采样点位的空间信息，是使用GPS记录采样点位的空间信息。所述的空间信息包括坐标信息和地面高程信息。

Claims (8)

1.一种表层岩溶带空间蓄水能力的估算方法，其特征在于，包括下述步骤：

a、利用探地雷达获得表层岩溶带雷达信息，同时记录采样点位的空间信息；

b、对所述的表层岩溶带雷达信息进行滤波处理，利用层位追踪的方法提取出表层岩溶带雷达信息中的深度信息；

c、根据空间信息和深度信息两者存在的空间体积差值获取表层岩溶带的体积；

d、获取表层岩溶带的裂隙率，计算出表层岩溶带的裂隙体积；核算表层岩溶带裂隙中的被土壤充填体积与未被土壤充填体积所占比例，计算被土壤充填裂隙的持水体积，所述的被土壤充填裂隙的持水体积与未被土壤充填体积之和即为表层岩溶带的蓄水体积；

e、步骤d所述的表层岩溶带的蓄水体积即为表层岩溶带空间蓄水能力。

2.根据权利要求1所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤a中，所述的利用探地雷达获得表层岩溶带雷达信息，包括下述步骤：

a1、设置数据采集窗口；

a2、通过探地雷达的发射天线与接收天线探测，在数据采集窗口反映出原始雷达影像；

a3、对比原始雷达影像所反映出的表层岩溶带的剖面信息，确定表层岩溶带的介电常数及层位划分原则，得到表层岩溶带深度信息。

3.根据权利要求2所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤b中，所述的进行滤波处理，是利用Reflex处理软件进行滤波处理，具体包括下述步骤：

b1、导入数据，设置表层岩溶带的介电常数；

b2、一维滤波/去直流漂移，除去零点漂移，设置两个参数，所述的两个参数分别为1time.ns和2time.ns；所述的1time.ns数值为时间窗口2/3处的值；所述的2time.ns数值为从最后样点倒数第5个样点处的数值；

滤波公式为：

其中，ω₁(t)是处理前第1道扫描回波与ω₁′(t)是处理后第1道扫描回波，M是扫描回波的总道数，ω_l(t)是处理前第l道扫描回波；

b3、静校正/移动开始时间，通过分析得到的第一个波峰处为直达波时间，把直达的开始位置移到零点处；

b4、增益/能量衰减，对探地雷达信号放大，选择所需增益放大的比例因子；

b5、二维滤波/抽取平均值

设置输入开始时间start time和结束时间end time；滤波公式为：

计算，其中Wm,n(t)为多次测量接收到的信号，m为各点测量序号，n为重复次数，M为扫描回波的道数；

b6、一维滤波/巴特沃斯带通滤波，设置低切频率与高切频率，得到有效波的图像；滤波公式为：

其中，n为滤波器阶数，ω_c为截止频率，ω_p为通频带边缘频率；H(ω)为振幅；

b7、二维滤波/滑动平均；滤波公式为：

其中n为滤波器阶数，ω_c为截止频率，ω_p为通频带边缘频率；H(ω)为振幅。

4.根据权利要求1或2所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤c中，所述的空间体积差值，是步骤a所述的空间信息中的地面高程信息与步骤b所述的深度信息存在的空间体积差值。

5.根据权利要求4所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤c中，所述的空间体积差值，按下述步骤计算：

c1、选取基准面作为计算表层岩溶带的体积的底界面，并在探地雷达获得表层岩溶带雷达信息后，运用空间插值，获得表层岩溶带的表层岩溶带下界；

c2、分别计算基准面与地面、基准面与表层岩溶带下界所构成的空间结构的体积，两者之差为空间体积差值，即表层岩溶带的体积。

6.根据权利要求5所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤c1中，所述的选取基准面，是选取低于表层岩溶带下界的基准面。

7.根据权利要求1或2所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤d中，所述的被土壤充填裂隙的持水体积为土壤田间持水度与土壤本身含水量的差值。

8.根据权利要求1或2所述的表层岩溶带空间蓄水能力的估算方法，其特征在于，步骤a中，所述的记录采样点位的空间信息，是使用GPS记录采样点位的空间信息。

Images