一种利用探地雷达获取田块尺度土体构型信息的勘测方法

1.本发明属于计量土壤学技术领域，具体涉及一种利用探地雷达获取田块尺度土体构型信息的勘测方法。


背景技术：

2.土体构型是各土壤层次的有序排列，是单个或多个不同时期的成土过程组合的结果。由于土体构型对土壤水分、养分以及污染物等运移具有决定性的作用，掌握土体构型空间分布信息对在农业、环境等领域具有重要意义。在实际应用中，针对不同尺度的土体结构信息，其勘测方式和精度要求不同。在区域尺度测绘和调查中，土体结构的调查通常是地学统计性为主。通过在调查区代表性点位进行破坏型采样获取一系列剖面信息，采样数据经处理后可以结合卫星数据以及其他相关观测数据，利用数据融合等方法构建研究相关预测模型，从而绘制区域土体构型分布图。而田块尺度土体构型勘测通常是为了获取其空间分布的精准信息，从而为精准灌溉、养分或污染物运移等实际应用提供支撑。此类勘测目的通过传统的高密度破坏性采样和空间插值显然是不能实现，尤其是空间变异较大土体构型分布情况。近二十年兴起的地球物理方法如探地雷达(gpr)等无损勘测手段，为获得高精度土体构型空间分布信息获取提供了可能。但目前相关的应用较少，且在勘测方法缺乏可靠的手段。
3.探地雷达探测土体构型的基本原理为：通过探地雷达向地下发射电磁波，其在土层界面处产生反射后被雷达接收，然后利用电磁波的传播时间进行反演计算获得土体层次信息。其中相关反演计算方法已相对成熟，而土体构型空间分布整体勘测效率较低且不确定性较大。通常利用探地雷达勘测土体构型空间分布的主要流程包括平行布线覆盖待测田块采样和空间插值重构两部分。由于探地雷达在测线方向上空间分辨率可高达厘米级，但测线间分辨率多在米级以上。在实际勘测中由于缺少必要的土体构型空间分布信息，平行测线的间距选择通常取决于调查者主观判断。然而若平行布线过密则工作量太大，若过于稀疏则采样的空间代表性有限。这会导致土体构型空间分布勘测结果的不确定性较大。此外，不合适的勘测方案会导致空间采样的均衡性交叉，从而降低常规空间插值方法的精度。
4.通过上述分析可知，探地雷达勘测方案对土体构型空间分布信息获取至关重要。原则上，空间变异大的方向需要相对较高的采样密度，而空间采样的均衡性直接决定着空间重构的精度。若既能充分利用探地雷达测线方向上高密度采样优势，又能选择合适的平行测线间距，则可以确保探地雷达采样在空间上的均衡性。有鉴于此，在利用探地雷达获取田块尺度高精度土体构型信息的勘测方法上，仍有进一步改进的空间。


技术实现要素：

5.针对现有技术中的上述不足，本发明提供一种利用探地雷达获取田块尺度土体构型信息的勘测方法。
6.为了达到上述发明目的，本发明采用的技术方案为：
7.提供一种利用探地雷达获取田块尺度土体构型信息的勘测方法，包括以下步骤：
8.s1、选择设备与配置参数：根据待勘测田块的土壤类型、土体结构及测试环境确定探地雷达系统的设备配置，根据数据需求完成探地雷达基础参数配置；
9.s2、采集初步勘测数据：选择3～5条中心交叉的测线，将探地雷达系统沿各条中心交叉的测线方向平移，完成对待勘测田块的初步初勘测，得到初步土体构型信息；
10.s3、确定详细勘测布线方向：利用初步勘测数据获取每条测线的微起伏指数mri，基于微起伏指数mri得到各条初步勘测测线对应的粗糙度，选择其中粗糙度最大的测线对应的方向作为详细勘测的布线方向；
11.s4、确定详细勘测布线间距：提取初步勘测数据中垂直于详细勘测布线方向上的土体构型，提取土体构型变异特征，获得待勘测田块土体构型的主体变异尺度序列，选择其中符合勘测需求的最大变异尺度作为正式勘测的关注尺度，将该关注尺度的1/4作为详细勘测中平行测线的布线间距；
12.s5、设计与实施详细勘测方案：根据详细勘测布线方向与间距制定出详细勘测布线方案，使用探地雷达系统沿按详细勘测方案进行勘测，计算每条测线上土体构型；
13.s6、重构土体构型三维空间：基于探地雷达勘测结果，对土体构型空间分布进行三维重构，提高数据分辨率，获得被调查田块高精度土体构型空间分布信息。
14.进一步地，步骤s2和s5中，在测线上应用共偏移距法co和宽角度反射折射法warr，将探地雷达系统沿各条中心交叉的测线方向平移，获取测线对应土体构型信息d，将获取到的土体构型使用savitzky-golay平滑滤波器进行低通滤波，获得测线上的降噪后土体构型信息d'。
15.进一步地，根据公式：
[0016][0017]
获取平均波速v；其中t为信号传播时间，t0为天线零间距时信号的理论传播时间，x为天线间距；
[0018]
根据公式：
[0019][0020]
获取土体构型信息d；其中t
co
为土体结构反射信号反射时间。
[0021]
进一步地，步骤s3中，根据公式：
[0022][0023]
获取微起伏指数mri；其中d'(x)与f(x)分别为降噪后的土体构型深度以及将其线性拟合后的对应结果，l为测线长度，f为地形波动频率。
[0024]
进一步地，步骤s4中，使用以morlet小波为母小波的小波分析方法获取主体变异尺度，具体方法包括以下子步骤：
[0025]
s4-1、根据公式：
[0026]
[0027]
获取小波系数wf；其中a与b分别为母小波(ψ)中尺度放缩与平移位置的控制系数；
[0028]
s4-2、根据公式：
[0029]
var(a)＝∫|wf(a，b)|2db
[0030]
获取小波方差var；
[0031]
s4-3、提取小波方差的极大值得到主体变异尺度。
[0032]
进一步地，步骤s6中，利用克里金插值法对土体构型空间分布进行三维重构；具体方法包括以下子步骤：
[0033]
s6-1、根据公式：
[0034][0035]
获取半方差函数γ；其中n为采样点的总数，h为采样点间距；
[0036]
s6-2、根据公式：
[0037][0038]
获取各已知点对应未知点的权重λ；
[0039]
s6-3、根据公式：
[0040][0041]
对未知点的值进行插值。
[0042]
本发明的有益效果为：
[0043]
1.本发明创新性地加入了初步勘测的步骤，通过少量先验信息获得土体构型空间变异性的主体性特征(方向和尺度)。基于此主体特征参数设置的详细勘测方案，使采样数据的空间分布与变异特征相适应，大大提高了采样点布局的合理性，一定程度上降低了土体构型三维空间重构对插值方法的依赖性。
[0044]
2.本发明首先确定详细勘测布线方向，进而确定合理的布线间距。本方法有效避免了现有常规勘测流程的盲目性，大大降低常规勘测由于采样代表性不足或过密而造成的土体构型三维空间重构精度的不确定性。
附图说明
[0045]
图1为本发明勘测方法流程图；
[0046]
图2为本发明实施例使用雷达系统相关warr与co原理示意图，其中tx与rx分别为
发射天线与接收天线；
[0047]
图3为本发明实施例中交叉测线pl1-4进行初步勘测的示意图，其中pl1-4的mri分别为1.84
×
10-3
、1.99
×
10-3
、1.73
×
10-3
与6.2
×
10-4
；
[0048]
图4为本发明实施例中垂直布线方向小波分析结果，其中选取的主体变异尺度为18.1m；
[0049]
图5为本发明实施例的布线方案与重构结果展示图；
[0050]
图6为传统布线方案与重构结果展示图。
具体实施方式
[0051]
下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
[0052]
实施例
[0053]
以某田块(36m
×
48m)土体构型勘测为例：
[0054]
参照图1，一种利用探地雷达获取田块尺度土体构型信息的勘测方法，包括以下步骤：
[0055]
步骤一：结合待勘测田块的土壤类型、土层厚度及土壤含水率等信息选择适宜的探地雷达天线频率；受雷达电磁波穿透深度和深度分辨率的影响，通常条件下对于1m以内土体构型勘测，应选用400mhz到900mhz的天线；本实施例中土体构型上部为沙壤土，土层厚度约1m，含水率约12％，因此选择了400mhz的探地雷达天线频率；由于本次勘测结果需要用于精细土壤水文过程研究，对于分辨率要求高，因此测线方向采样分辨率设置为0.05m。
[0056]
步骤二：从田块中选择三个代表性点位，使用宽角反射折射法确定土层中的平均波速；进一步地，选择4条中心交叉的测线(pl1-4)作为初步勘测布线方案(图3)；使用共偏移距法分别沿四条中心交叉的测线方向进行勘测，获取各个方向土层厚度分布信息(图2)。
[0057]
优选的，将获取到的土体构型使用savitzky-golay平滑滤波器进行低通滤波，获得测线上的降噪后土体构型信息d'。
[0058]
具体地说，根据公式：
[0059][0060]
获取平均波速v；其中t为信号传播时间，t0为天线零间距时信号的理论传播时间，x为天线间距；
[0061]
根据公式：
[0062][0063]
获取土体构型信息d；其中t
co
为土体结构反射信号反射时间。
[0064]
步骤三：利用步骤二中的初步勘测结果计算每条测线的微起伏指数(mri)，得到各条初步勘测测线对应的粗糙度(图3)。pl1-4分别为1.84
×
10-3
、1.99
×
10-3
、1.73
×
10-3
与6.2
×
10-4
；将其中粗糙度最大的测线(即pl2)对应的方向设置为详细勘测方案的布线方向，
即正式勘测布线方向应设置为60
°
。
[0065]
具体地说，根据公式：
[0066][0067]
获取微起伏指数mri；其中d'(x)与f(x)分别为降噪后的土体构型深度以及将其线性拟合后的对应结果，l为测线长度，f为地形波动频率。
[0068]
找到其中粗糙度最大的初步勘测测线对应的方向并将该方向设置为详细勘测的布线方向，在本实例中选择pl2对应方向。
[0069]
步骤四：选择垂直于详细勘测方向(pl2)的测线数据(即示例中pl4)作进一步的土体构型空间变异特征分析；利用小波分析方法提取该勘测方向上土体构型的主要变异特征数据(如图4中各峰值对应的变异尺度)；根据勘测精度的要求选择关注的主体变异尺度，并将该关注尺度的1/4作为详细勘测方案中平行测线的布线间距(图5)，本实施例中关注的主体变异尺度为18.1m，布线间距设置为4m。
[0070]
具体地说，使用以morlet小波为母小波的小波分析方法获取主体变异尺度，具体方法包括以下子步骤：
[0071]
根据公式：
[0072][0073]
获取小波系数wf；其中a与b分别为母小波(ψ)中尺度放缩与平移位置的控制系数；
[0074]
根据公式：
[0075]
var(a)＝∫|wf(a，b)|2db
[0076]
获取小波方差var；
[0077]
提取小波方差的极大值得到主体变异尺度。
[0078]
步骤五：根据步骤三和步骤四提供的详细勘测布线方向与间距制定出详细勘测布线方案，即沿pl2方向布置13条间距为4m的平行测线，使用共偏移距法，沿侧线进行探地雷达勘测作业并计算每条测线的土层厚度分布数据。
[0079]
步骤六：利用普通克里金插值法对步骤五中获得土体构型数据进行三维空间重构，最后获得被调查田块高精度土体构型空间分布信息，如图5的分辨率可达5cm。
[0080]
具体地说，根据公式：
[0081][0082]
获取半方差函数γ；其中n为采样点的总数，h为采样点间距；
[0083]
根据公式：
[0084][0085]
获取各已知点对应未知点的权重λ；
[0086]
根据公式：
[0087][0088]
对未知点的值进行插值。最后获得被调查田块高精度土体构型空间分布信息(图5)。
[0089]
本发明可以相对传统勘测方案，以布线方向为150
°
、布线间距4m为例(图6)，在相同工作量(即布线间距相同)的情况下，本勘测方法有效规避了盲目性高的弊端；在相同的布线间距下，重构质量大大提高；这也说明，在无初步勘测的指导下，如想达到相同的三维重构质量，需要更小的布线间距。
[0090]
本发明创新性地加入了初步勘测的步骤，通过少量先验信息获得土体构型空间变异性的主体性特征(方向和尺度)。基于此主体特征参数设置的详细勘测方案，使采样数据的空间分布与变异特征相适应，大大提高了采样点布局的合理性，一定程度上降低了土体构型三维空间重构对插值方法的依赖性。有效提高土层厚度勘测工作效率。
[0091]
本方法有效避免了现有常规勘测流程的盲目性，大大降低常规勘测由于采样代表性不足或过密而造成的土体构型三维空间重构精度的不确定性。
[0092]
于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
[0093]
此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。