一种基于探地雷达的红壤关键带地下结构空间预测方法与流程
本发明涉及一种基于探地雷达的红壤关键带地下结构空间预测方法，属于计量土壤学技术领域。
背景技术：
地球关键带是地球表层系统中土壤圈与大气圈、生物圈、水圈、岩石圈物质迁移和能量交汇的区域，在垂直方向上包括从植被冠层到地下含水层的区域，是人类生存与发展的核心区域。土壤的形成经历了基岩的风化过程：未风化→微风化→中等风化→强风化→残积土。广义上，未风化的基岩层定义为关键带的下层边界，关键带的厚度也即从地表到未风化基岩的深度。该厚度对于区域甚至全球的农业、生态、环境保护、气候变化意义重大，例如区域储水能力、碳储量。将关键带的厚度、形状与地球物理过程相结合是关键带科学重要突破口之一，理解关键带各结构层次的厚度及其空间分布特征是气候变化、生态管护、水资源安全、自然灾害防治工作中日趋活跃的核心技术难题。
红壤是我国南方典型的土壤类型，分布面积约为220万平方公里，约占国土面积的23％，支撑了我国40％的人口。雨热条件优越，经济发展较快，农田集约化程度高，农业产值占全国的1/2。当前高强度的人为利用和日益增加的酸沉降使得该地区水土流失、土壤酸化以及养分贫瘠等土壤退化问题极其严重，导致该地区生态与环境遭受严重破坏。土壤安全是土壤可持续发展的系统战略框架。因此，红壤关键带与区域地理尺度、景观特征是密切相关的，通过对红壤关键带不同结构层次的反演及空间分布的预测，能够为跨学科的关键带工作提供基础的数据支撑与技术支持，进而为生态过程、生物地球化学过程和水文过程等方面相关工作的推进提供空间参考。
红壤区地下结构的空间变异既有结构性因素，又有随机性因素，结构性因素变异的空间尺度较大。关键带是具有高度空间异质性的非均质体，其结构与理化性质垂直分布特征受地形、地质、植被、气候等多因素的综合影响，这种垂直分布规律不仅为不同景观区域内土壤碳、氮等元素储量合理评测提供了理论依据，也为三维关键带模拟提供了理论支持。数字土壤制图理论与方法可为关键带三维空间分布特征表达技术提供重要借鉴。空间变异现象可以采用定量分类、多元统计、模糊分类、分形理论、数量地貌学等方法来描述。数字土壤制图方法按照分析性质可以分为基于数理统计方法，如线性模型、地理加权回归等；基于地统计方法，如普通克里格、协同克里格等；基于模糊逻辑方法，如模糊c均值聚类、模糊系统推理等；基于神经网络方法，如bp网络、径向基神经网络等；基于机器学习方法，如决策树、支持向量机等；这些算法在不同条件下预测性能相差较大。增强环境变量的特征表达能力并消除冗余信息是土壤制图领域另一技术热点。三维土壤制图一直是计量土壤学的重要主题。国内外学者已针对三维土壤制图提出了众多可行的技术方案，如3d克里金、神经网络、回归树等。
在地下关键带结构反演方面，现有地球物理勘探技术能够提供一定的支撑。探地雷达是由发射天线向地下介质发射电磁脉冲波，电磁脉冲波在地下介质中传播时，遇到介质中的电磁性(电阻率、介电常数及磁导率)差异界面会发生反射、透射等。被反射的电磁波传回地表，由接收天线接收。使用低频探地雷达能够分析雷达波的频谱特征，并反演地表至未风化基岩层次间不同风化程度。不同风化程度岩石的物理性质描述方面已具备大量指标，能够准确模拟其颜色、岩性、结构、构造、风化程度的差异，如天然密度、天然含水量、土粒密度、天然干密度、孔隙比、最大干密度、最佳含水量、孔隙率、饱和度等。三维地质建模现已取得重要进展，有限的地质剖面数据与预测模型的有机结合是主要技术手段。
然而，常规的探地雷达工作仍是基于测线的方式，工作效率较为低下，无法预测生成大区域地区的地下关键带结构。另外，地球物理勘探的结果在三维模型表达方式、空间异质性表达等方面也无法与应用部门的相关技术无缝衔接。有鉴于此，红壤关键带地下结构空间预测技术迫切需要解决地球物理勘探与基于关键带-景观的空间预测方法的有机结合问题，仍需解决一系列技术问题：
(1)现有工作的地理尺度需要进行外推。目前，全球关键带监测网络已建成64个节点。通过大量的原位测定实验，关键带科学观测与模拟主要集中在田间尺度。例如，我国在江西省鹰潭市建立的红壤关键带观测流域面积仅为50公顷。通过大量的地质钻井，已获取了该流域红壤关键带地下结构的相关信息，但是市级、省级及其他地区的红壤区的地下结构信息仍然较为匮乏。
(2)红壤关键带结构、厚度定义尚不够明确，缺乏统一的标准。地质学仅使用风化系数与波速比等指标定量化岩石风化程度；土壤学重点关注成土母质以上生长植物的疏松表层；水文地质学主要研究地下水的数量和质量及其时空变异规律。因此，需要将关键带作为整体框架，从地质学、土壤学、水文学、生态学等不同领域学科综合视角出发，基于对关键带复杂物理过程、化学过程、生物过程的深入剖析，厘定关键带厚度划分标准。
(3)计算机理问题。土壤制图方法多为知识发现模型与土壤发生学简单结合后的案例应用，相对忽视了关键带各结构层次间的能量转换与物质迁移，缺乏系统的方法创新与集成技术。常规的土壤调查采样深度为1米左右，红壤关键带厚度约为10-30米。钻井的高成本导致获取的数据非常有限，且关键带结构的空间变异较高，导致预测模型的不确定性较高。现有三维土壤制图方法较多关注一米深度土壤预测，较少关注关键带不同风化层理化属性垂直维变异特征。
(4)现有计量土壤学技术较多关注二维制图或植物根系生长层为主的浅表土体(1-2m)三维制图，缺乏结合地球物理勘探方法的关键带结构、厚度反演，尚未综合、系统地开发表达关键带空间结构的相关技术体系，尤其是面向地下关键带的不规则三维空间模拟与表达。
上述技术缺陷已严重影响了红壤区地下结构信息的高效获取，更是直接阻碍了不同行业部门的相关工作，降低了生产加工土壤信息产品的效率与精度，也间接给国家经济规划与调控带来了潜在的损失。
技术实现要素：
本发明所要解决的技术问题是提供一种基于探地雷达的红壤关键带地下结构空间预测方法，不仅能够给出详细的关键带结构定义、数据库构建、探地雷达作业流程，同时也能定量化关键带各层次结构厚度/类型与环境变量的空间耦合关系。
本发明为了解决上述技术问题采用以下技术方案：本发明设计了一种基于探地雷达的红壤关键带地下结构空间预测方法，包括如下步骤：
步骤a.根据预设数据结构类型，构建目标红壤关键带地下结构数据库；
步骤b.采用探地雷达，沿预设测线，针对目标红壤关键带进行测量，获得目标红壤关键带雷达影像数据；
步骤c.根据目标红壤关键带地下结构数据库，对目标红壤关键带雷达影像进行解译，获得雷达影像解译结果；
步骤d.基于雷达影像解译结果，更新目标红壤关键带地下结构数据库；
步骤e.基于目标红壤关键带地下结构数据库，结合环境变量数据库，构建目标红壤关键带地下结构三维预测模型，进而获得目标红壤关键带地下不同层次结构空间分布。
作为本发明的一种优选技术方案，所述步骤a包括如下步骤：
步骤a1.根据目标红壤关键带的历史土壤数据，获得目标红壤关键带中各个历史样点的坐标信息，以及历史样点土壤剖面的分层信息；其中，分层信息包括各层土壤的类型信息、深度信息、厚度信息；
步骤a2.针对目标红壤关键带布设野外钻井样点，并获得各个野外钻井样点的坐标信息；
步骤a3.分别针对各个野外钻井样点进行钻探，获得各野外钻井样点位置、自地表至基岩之间土壤剖面的各分层信息；
步骤a4.针对目标红壤关键带中各个历史样点、各个野外钻井样点的坐标信息，以及土壤剖面的分层信息，基于预设统一坐标系，按预设数据结构类型，进行数据标准化，进而构建目标红壤关键带地下结构数据库。
作为本发明的一种优选技术方案：所述步骤a2中，根据目标红壤关键带的区域大小、资金情况、土地利用面积与道路可达性，使用目的性采样方法，针对目标红壤关键带布设野外钻井样点。
作为本发明的一种优选技术方案，所述步骤c包括如下步骤：
步骤c1.针对目标红壤关键带雷达影像，进行静校正/移动开始时间操作，消除探地雷达探地测量操作过程中未测量地表以下部分的时间，更新目标红壤关键带雷达影像；
步骤c2.针对目标红壤关键带雷达影像，进行去直流漂移操作，消除零点漂移，更新目标红壤关键带雷达影像；
步骤c3.针对目标红壤关键带雷达影像，滤除其中的干扰噪声数据，更新目标红壤关键带雷达影像；
步骤c4.针对目标红壤关键带雷达影像，进行背景去除操作，消除影像中的水平部分，更新目标红壤关键带雷达影像；
步骤c5.针对目标红壤关键带雷达影像，进行雷达信号增强处理，将预设深度标准的雷达影像信号进行放大；
步骤c6.根据目标红壤关键带地下结构数据库中各个历史样点、各个野外钻井样点的坐标信息，以及土壤剖面的分层信息，针对目标红壤关键带雷达影像进行对比匹配，获得雷达影像解译结果。
作为本发明的一种优选技术方案：所述步骤c3中，采用带通滤波、高通滤波、低通滤波中的任意一种方法，针对目标红壤关键带雷达影像进行滤波处理，滤除其中的干扰噪声数据，更新目标红壤关键带雷达影像。
作为本发明的一种优选技术方案：所述步骤c5中，采用反射回波幅度变换技术和振幅处理技术，针对目标红壤关键带雷达影像，进行雷达信号增强处理，将预设深度标准的雷达影像信号进行放大。
作为本发明的一种优选技术方案，所述步骤d包括如下步骤：
步骤d1.根据预设预测结果空间分辨率，确定从雷达影像解译结果中反演探测样点之间的水平间距；
步骤d2.针对雷达影像解译结果，将沿测线的连续测量数据，转换为以水平间距为间隔的各个反演探测样点；
步骤d3.根据水平间距，获得各个反演探测样点的坐标信息，并基于雷达影像解译结果，分别获得各个反演探测样点土壤剖面的分层信息；
步骤d4.根据雷达影像解译结果所对应各个反演探测样点土壤剖面的分层信息，更新目标红壤关键带地下结构数据库。
作为本发明的一种优选技术方案，所述步骤d1包括如下步骤：
步骤d1-1.根据预设预测结果空间分辨率res，按如下公式获得：
获得雷达影像解译结果所对应测线上反演探测样点的数目n_total，其中，m表示测线的长度，n表示测线上土壤剖面的层数；
步骤d1-2.按m/n_total，确定从雷达影像解译结果中反演探测样点之间的水平间距。
作为本发明的一种优选技术方案，所述步骤d3中，基于雷达影像解译结果，分别针对各个反演探测样点，获得反演探测样点，以及该反演探测样点为中心、结合前后预设数目其它反演探测样点对应土壤剖面分层信息的平均值，作为该反演探测样点土壤剖面的分层信息。
作为本发明的一种优选技术方案，所述步骤e包括如下步骤：
步骤e1.针对目标红壤关键带地下结构数据库中的各个样点，获得样点分别对应预设各类环境变量的信息；
步骤e2.分别针对土壤剖面的各个土壤层，以各个样点分别所对应预设各类环境变量信息为自变量，以土壤剖面土壤层的厚度为因变量，采用逐步回归方法，筛选获得土壤剖面土壤层的厚度所对应的最优环境变量；
步骤e3.分别针对土壤剖面的各个土壤层，以土壤层的厚度、以及该土壤层所对应的最优环境变量作为输入数据，针对预设预测模型进行训练，获得目标红壤关键带地下结构各土壤层厚度预测模型；
同时，以各个样点分别对应预设各类环境变量的信息为自变量，各个样点分别所对应土壤剖面的总深度为因变量，针对预设预测模型进行训练，获得目标红壤关键带地下结构土壤剖面深度预测模型；
步骤e4.采用目标红壤关键带地下结构各土壤层厚度预测模型进行预测，获得目标红壤关键带地下结构各土壤层厚度分布；同时，采用目标红壤关键带地下结构土壤剖面深度预测模型进行预测，获得目标红壤关键带地下结构土壤剖面深度分布，并构建土壤剖面深度分布图层；
步骤e5.将目标红壤关键带地下结构各土壤层厚度分布，以及土壤剖面深度分布图层，保存为栅格文件，获得目标红壤关键带地下不同层次结构空间分布。
本发明所述一种基于探地雷达的红壤关键带地下结构空间预测方法，采用以上技术方案与现有技术相比，具有以下技术效果：
(1)本发明所设计一种基于探地雷达的红壤关键带地下结构空间预测方法，涉及到生态过程、水文过程、生物地球化学过程及其相互耦合关系，创新性地通过风化层的地球物理勘探与关键带科学研究的有机结合，从岩石风化--土壤形成这一土壤发生学与地球化学视角，解决了红壤区关键带地下结构划分与数据库构建的一系列技术难题，为相关科学研究与工程作业提供了技术参考；
(2)本发明所设计一种基于探地雷达的红壤关键带地下结构空间预测方法，具有较好的迁移性，其中，地下结构层次的空间预测技术是区域级关键带结构制图的关键。本发明提出的预测方法具有“高内聚、低耦合”的技术特点，部分的技术细节能够较为方便地移植至其他相关的技术研发，根据不同数据集开发的关键带地下结构预测机制，能够为面向地球表层系统科学及相关领域的三维制图，提供更广泛的知识迁移能力；
(3)本发明所设计一种基于探地雷达的红壤关键带地下结构空间预测方法，采用关键带地下结构与环境变量的定量化机制，三维关键带制图技术障碍之一在于如何量化不同结构层次与环境变量的关系，空间预测方法也受限于各层次垂直维/水平维的空间变异，定量化大区域关键带结构与环境变量的耦合关系也易受到空间尺度、误差传播等因素的制约，基于此，本发明重点关注了充分利用现有历史土壤数据、土壤专家知识、高维环境变量，并建立了一整套关键带层次结构的空间预测流程，为水资源安全、生态系统管护、气候变化等重大问题的解决提供了基础数据支持。
附图说明
图1是本发明所设计基于探地雷达的红壤关键带地下结构空间预测方法的流程图；
图2是探地雷达共中心点法测量示意图；
图3是探地雷达数据解译的步骤流程图；
图4是将连续的层位深度数据转换为固定间隔的点位数据示意图；
图5是基于环境变量数据库构建三维预测模型的步骤流程图；
图6是典型环境变量与结果分析流程的步骤流程图；
图7是本发明实施例中探地雷达影像解译结果；
图8是本发明实施例中均质红土层厚度的空间分布图；
图9是本发明实施例中网纹红土层厚度的空间分布图；
图10是本发明实施例中风化层厚度的空间分布图；
图11是本发明实施例中总厚度预测的空间分布图；
图12是本发明实施例中总厚度预测结果减去所有层次预测结果累加和结果的空间分布图；
图13是本发明实施例中典型序列的空间位置；
图14是本发明实施例中典型序列不同结构层次厚度的空间分布图。
具体实施方式
下面结合说明书附图对本发明的具体实施方式作进一步详细的说明。
如图1所示，本发明设计了一种基于探地雷达的红壤关键带地下结构空间预测方法，应用于实际当中，具体包括如下步骤：
步骤a.根据预设数据结构类型，构建目标红壤关键带地下结构数据库，然后进入步骤b。
我国南方第四纪红土的成因类型包括沉积型和风化残积型。前者根据沉积相的不同可分为坡积、冲积、风成等类型。残积型红土是南方红土区最常见的红土类型，主要包括花岗岩、红色砂页岩、石灰岩等基岩经历湿热风化后残留在原地的红土。因此，由地表至基岩，红壤关键带的结构也相差较大，部分典型的红土剖面地层结构自地表至基岩如：均质红土→网纹红土→风化层→基岩、黄棕色土→网纹黄棕色土→网纹红土→风化层→基岩、黄棕色土→网纹红土→风化层→基岩。
如果红壤区的岩石类型不同或地形起伏较为复杂，红土剖面的层次差异性更大。因此，为了便于预测模型格式化输入参数，使用如下结构对各层次进行表达：layer＝{layer_1,layer_2,…,layer_i,…,layer_i}，其中layer_i表示第i层的类型信息、深度信息、厚度信息，i表示该红土剖面的层次数。
注：弱网纹发育层质地黏重、较硬，最大的特点是网纹发育不明显。考虑到弱网纹发育层与网纹红土层质地、结构差异性不够显著，将弱网纹发育层的厚度合并至网纹红土层。
上述步骤a应用中，具体包括如下步骤：
步骤a1.根据目标红壤关键带的历史土壤数据，获得目标红壤关键带中各个历史样点的坐标信息，以及历史样点土壤剖面的分层信息；其中，分层信息包括各层土壤的类型信息、深度信息、厚度信息。
实际应用中，对于目标红壤关键带中的各个历史样点，使用剖面位置的市县名称的首字母组合，给每一个历史样点设定唯一的id，如果该市县具有多个剖面，需要在字母后面添加数字，数字的位数需要根据剖面数量的位数决定。
示例：南昌市具有362个样点，则命名规则是“nc_***”，该地区的红土剖面编号如下：nc_001，nc_002，…，nc_362。
步骤a2.根据目标红壤关键带的区域大小、资金情况、土地利用面积与道路可达性，使用目的性采样方法，针对目标红壤关键带布设野外钻井样点，并获得各个野外钻井样点的坐标信息；其中，野外钻探的之前，需要使用手持gps记录钻井的经纬度信息，即获得各个野外钻井样点的坐标信息，钻井作业的具体作业标准参见中华人民共和国地质矿产行业标准-地质岩心钻探规程(dz/t0227-2010)。
步骤a3.分别针对各个野外钻井样点进行钻探，获得各野外钻井样点位置、自地表至基岩之间土壤剖面的各分层信息，具体来说，在钻探的过程中，样品会间歇地从钻头中取出，需要有一名有经验的地质学工作者观察每一次取出来的样品，并立刻记录不同层次结构的出现深度。
步骤a4.针对目标红壤关键带中各个历史样点、各个野外钻井样点的坐标信息，以及土壤剖面的分层信息，基于预设统一坐标系，统一转换成地理坐标(gcs_wgs_1984)，并与投影坐标(wgs_1984_albers)，将每一个样点的结构的厚度转换标准的国际制(m)，按预设数据结构类型，进行数据标准化，进而构建目标红壤关键带地下结构数据库。
针对目标红壤关键带中所有历史样点、所有野外钻井样点，统一划归为样点，针对所有样点，选取20％作为验证样点集合test_obser，其余的样点为训练样点集合，记为train_obser。
步骤b.采用探地雷达，沿预设测线，针对目标红壤关键带进行测量，获得目标红壤关键带雷达影像数据，然后进入步骤c。
上述步骤b中，具体执行如下步骤：
步骤b1.野外测线布置。野外测线布置需要制定勘测方案，并选择合适的低频天线。红壤关键带地下结构包含了若干地质层位，因此探地雷达的探测对象特点显著。根据野外钻井的信息描述，需要根据基岩的埋深来决定天线的探测深度。一般情况下，红壤关键带地下结构的水平起伏形状与地表形状呈平行状。测线布设需要重点关注探测目标水平尺度的大小及水平分辨率。测线的间距应该同时小于或等于目标尺度与分辨率尺度。
测线的选择应根据调查的工作时限、土地利用条件、天气状况综合决定。建议在土地利用条件允许的情况下，尽可能布置较长的测线，以便获取更多的地下结构信息。
步骤b2.探地雷达发射、接收天线频率选择。天线是雷达发射与接收电磁波的部件。天线中随时间变化的电荷或电流在天线的周围激发电磁波，电磁波可以脱离天线传向远处，形成电磁辐射。探地雷达的频率可大致分为低频、中频与高频三种类型。记电磁波的波长为λ，各类型的雷达天线特征如下：
低频：80mhz以下的为低频天线，通常采用非屏蔽式半波偶极子杆状天线。天线的总长度等于λ/2，发射频率低，介质中衰减小，能够探测埋藏较深的目标。
中频：频率在100mhz-1ghz范围内的天线，采用屏蔽式半波偶极子天线，天线总长度为λ/2。
高频：频率高于1ghz，有1.0ghz、1.5ghz、2.0ghz多种。路面检测使用的高频天线常采用喇叭形状，该天线辐射能量集中，分辨率高，以提高辐射效率。工程检测的高频天线仍采用屏蔽半波偶极子天线，体积小，便于操作。部分频率的雷达探测深度与应用领域见下表1。
表1
步骤b3.野外探地雷达作业。根据作业场地实际情况，可以采用“离散”测量(点测)与“连续”测量(线测)两种工作方式。点测适宜于地表起伏较为剧烈的情况，发射天线与接收天线每隔一定的距离探测一次。线测适宜于地表起伏较为平缓的情况，连续地对每一个位置进行探测。
现场测量时需要对雷达参数进行初始化设置，具体设定的参数包括天线距、增益点数、时间窗口大小、介电常数大小、扫描速率、每秒扫描数等内容。其中，时窗设置涉及到雷达信号的采集深度，不仅要保证所需要的信息都在探测深度内，也不能设置过大而影响到垂直分辨率。选取探测深度d为目标深度的1.5倍。根据探测深度depth、波速v(或介电常数ε)确定采样时窗长度：
式中，width为采样时窗(ns)，depth为拟探测深度(m)，v是平均雷达波速(m/ns)。由该公式可以得出：时窗越大，探测深度也越大，分辨率也越低。
野外作业还需要使用gps仪器记录测线的位置信息。如果测线是条直线，则只需要记录起始位置信息即可。如果测线在作业的过程中发生了偏移或拐弯，则需要记录野外作业时测线特征点的位置信息(拐点)。同时，测量的过程中，还需要记录好土地利用、天气、土壤水分状况等信息，以便室内作业更加准确地进行雷达影像解译。
步骤b4.根据设备的情况，选择是否进行共中心点法(cmp法)测量，如果设备发射天线与接收天线是分离的，或具有两个集成式天线(发射天线与接收天线绑定在一个设备上)，则进行cmp测试：
在钻井周围地区选择具有代表性的区域，发射天线与接收天线向相反的方向移动，如图2所示，该技术借鉴了地震反射勘探中的多次覆盖技术，能够获取该样点的相对介电常数。
否则，进入步骤b5。
步骤b5.根据相对介电常数计算雷达波速(v)：
式中，ε为相对介电常数，c是电磁波在大气中的传播速度，约为3×108m/s。介电常数最大的物质是水(介电常数为81)，最小的是空气与金属(介电常数为1)。因为水是自然界中介电常数最大、电磁波速最低的介质，一般情况下，野外作业中地下结构介电常数的主要影响因素是水分含量。
如果已经执行cmp探测，则可以计算不同深度的介电常数，进而由公式(2)获取到雷达波速：
式中，tx是从发射到接收经过目标的反射波的双程行进时间，t0是零偏移距时的反射波双程行进时间，dis是发射天线和接收天线间的距离。
否则，可以参考不同介质类型来获取大致的相对介电常数，如下表2所示，进而结合钻井信息校正探测的深度信息。
表2
雷达波速的计算也可以通过其的方法计算，如已知目标深度法、层状反射体法、点源反射体法等。
步骤c.根据目标红壤关键带地下结构数据库，对目标红壤关键带雷达影像进行解译，获得雷达影像解译结果，然后进入步骤d。
上述步骤c实际应用中，如图3所示，具体包括如下步骤：
步骤c1.针对目标红壤关键带雷达影像，进行静校正/移动开始时间操作，消除探地雷达探地测量操作过程中未测量地表以下部分的时间，更新目标红壤关键带雷达影像。
步骤c2.针对目标红壤关键带雷达影像，进行去直流漂移操作，消除零点漂移，更新目标红壤关键带雷达影像。该技术是地震信号处理中较为常用的界面归位处理技术之一。由于探地雷达的反射波主要受到不同介质的典型差异影响，因此偏离目标测试位置的反射信息也都记录至接收终端。经过偏移处理的雷达剖面，更能真实反映地下介质的位置。
步骤c3.采用带通滤波、高通滤波、低通滤波中的任意一种方法，针对目标红壤关键带雷达影像进行滤波处理，滤除其中的干扰噪声数据，更新目标红壤关键带雷达影像。
步骤c4.针对目标红壤关键带雷达影像，进行背景去除操作，消除影像中的水平部分，更新目标红壤关键带雷达影像。
步骤c5.采用反射回波幅度变换技术和振幅处理技术，针对目标红壤关键带雷达影像，进行雷达信号增强处理，将预设深度标准的雷达影像信号进行放大。由于受到不同介质的干扰，雷达信号随着深度的增加出现一定程度的衰减，该步骤的主要目的是对雷达信号进行增强处理,改善图像质量以便进行识别。
步骤c6.根据目标红壤关键带地下结构数据库中各个历史样点、各个野外钻井样点的坐标信息，以及土壤剖面的分层信息，针对目标红壤关键带雷达影像进行对比匹配，获得雷达影像解译结果。
步骤d.基于雷达影像解译结果，更新目标红壤关键带地下结构数据库，然后进入步骤e。
实际应用中，上述步骤d包括如下步骤：
步骤d1.根据预设预测结果空间分辨率，确定从雷达影像解译结果中反演探测样点之间的水平间距。
其中，步骤d1具体包括如下步骤：
从雷达影像解译结果中，将连续的层位深度数据转换为固定间隔的点位数据，水平间隔为1m，生成数据文件：layer_i_depth＝{d0,d1,d2,…,dq,…,dq}，其中dq表示测线第qm处的样点，q表示测线距离，示意图如图4所示。
步骤d1-1.根据预设预测结果空间分辨率res，按如下公式获得：
获得雷达影像解译结果所对应测线上反演探测样点的数目n_total，其中，m表示测线的长度，n表示测线上土壤剖面的层数。
步骤d1-2.按m/n_total，确定从雷达影像解译结果中反演探测样点之间的水平间距。
步骤d2.针对雷达影像解译结果，将沿测线的连续测量数据，转换为以水平间距为间隔的各个反演探测样点。
步骤d3.根据水平间距，获得各个反演探测样点的坐标信息，并基于雷达影像解译结果，分别针对各个反演探测样点，获得反演探测样点，以及该反演探测样点为中心、结合前后预设数目其它反演探测样点对应土壤剖面分层信息的平均值，作为该反演探测样点土壤剖面的分层信息。
其中，根据水平距离，计算各个反演探测样点的经纬度信息。为了有效降低探地雷达影像的解译误差，将各个反演探测样点临近若干样点(n_neigh)的平均值赋给该反演探测样点的深度数值：
如果res为偶数，则di＝average(di-res/2，…，di-2，di-1，di，di+1，di+2，…，di+res/2)，di+1表示di反演探测样点沿测线方向移动1m后的反演探测样点，average()表示取平均值函数，di-res/2表示di样点朝测线起始方向移动res/2m的反演探测样点，反之，di+res/2表示di反演探测样点朝测线移动方向移动res/2m的反演探测样点。
如果res为奇数，则
如图4所示，如果res＝3，则衍生样点6个{d0_aver,d1_aver,d2_aver,d3_aver,d4_aver,d5_aver}，其中d0_aver＝average(d0,d1,d2)，d1_aver＝average(d1,d2,d3,d4,d5)。
计算每一个反演样点不同层位的厚度信息。至此步骤，所有的反演探测样点记录的均为层位的深度信息，记反演探测样点j的各层位的深度信息为depth_j＝{depth_j_1,depth_j_2,…,depth_j_k,…,depth_j_k}，k表示该红土剖面共k个层次，则第k层的厚度为：
thickness_j_k＝depth_j_k-depth_j_k-1(5)
式中，depth_i_k-1表示反演探测样点j的第k-1个层次的深度。
因此，k个层位的厚度为：thickness_j＝{depth_j_1,depth_j_2-depth_j_1,depth_j_k-depth_j_k-1,…,depth_j_k-depth_j_k-1}。式中，depth_j_k-depth_j_k-1表示第k层的深度减去第k-1层的深度。
针对反演探测样点，随机选取样点的20％作为验证样点集合test_new，其余80％的样点为模型训练样点，记为train_new。
步骤d4.根据雷达影像解译结果所对应各个反演探测样点土壤剖面的分层信息，更新目标红壤关键带地下结构数据库。
步骤e.基于目标红壤关键带地下结构数据库，结合环境变量数据库，构建目标红壤关键带地下结构三维预测模型，进而获得目标红壤关键带地下不同层次结构空间分布。
步骤e在实际应用中，如图5所示，具体包括如下步骤，典型的环境变量与结果分析流程如图6所示。
步骤e1.针对目标红壤关键带地下结构数据库中的各个样点，获得样点分别对应预设各类环境变量的信息。
根据关键带-景观模型定义，对于红壤关键带地下结构影响最为显著的五大要素分别是：地形、母质、生物、气候、时间。围绕上述影响因素，收集相应的环境变量作为预测不同层次厚度的自变量，例如数字高程模型、地质图、遥感图像、年均气候、年均降雨、土地利用与土壤类型等数据，典型的环境变量如图6所示。部分的数据，如地质图、土地利用与土壤类型图为矢量图层，需要转换成栅格图层，栅格图层的空间分辨率为res米。
步骤e2.分别针对土壤剖面的各个土壤层，以各个样点分别所对应预设各类环境变量信息为自变量，以土壤剖面土壤层的厚度为因变量(train_obser+train_new)，采用逐步回归方法，筛选获得土壤剖面土壤层的厚度所对应的最优环境变量。
步骤e3.分别针对土壤剖面的各个土壤层，以土壤层的厚度、以及该土壤层所对应的最优环境变量(train_obser+train_new)作为输入数据，针对预设预测模型进行训练，如图6所示，获得目标红壤关键带地下结构各土壤层厚度预测模型。其中，预设预测模型可以选用传统的数理统计方法、基于机器学习方法、基于地统计方法、基于神经网络方法等。
同时，以各个样点分别对应预设各类环境变量的信息为自变量，各个样点分别所对应土壤剖面的总深度为因变量，针对预设预测模型进行训练，获得目标红壤关键带地下结构土壤剖面深度预测模型。
步骤e4.采用目标红壤关键带地下结构各土壤层厚度预测模型进行预测，获得目标红壤关键带地下结构各土壤层厚度分布；同时，采用目标红壤关键带地下结构土壤剖面深度预测模型进行预测，获得目标红壤关键带地下结构土壤剖面深度分布，并构建土壤剖面深度分布图层。
并且可以根据验证样点集合的观测值与雷达影像解译结果对预测结果进行交叉验证，精度评价指标可选用皮尔逊相关系数、确定系数、均方根误差、平均误差等常规指标。
步骤e5.将目标红壤关键带地下结构各土壤层厚度分布，以及土壤剖面深度分布图层，保存为栅格文件，获得目标红壤关键带地下不同层次结构空间分布。
下面以江西省鹰潭市孙家流域的第四纪红土层厚度空间预测为例
孙家流域位于江西省鹰潭市，面积50.5ha，海拔高度在34～55m之间。流域内的土壤母质为第四纪红粘土与红砂岩，土壤类型为红壤。依据美国土壤系统分类方法，土壤类型为老成土(旱地土壤)和人为土(水稻田土壤)。该流域主要土地利用方式花生(48.7％)、稻田(24.8％)、橘园和葡萄园(19.8％)等。该区域是中国典型红壤关键带观测站，地下结构层次自上而下包括均质红土层、网纹红土层、风化层。基于低频探地雷达(60mhz)进行第四纪红土层厚度空间预测，具体步骤包括：
步骤a.根据预设数据结构类型，构建目标红壤关键带地下结构数据库，其中，红壤关键带地下结构定义：均质红土→网纹红土→风化层→基岩。
步骤b.采用探地雷达，沿预设测线，针对目标红壤关键带进行测量，获得目标红壤关键带雷达影像数据；其中，布置的测线距离共5km，探地雷达选择60mhz的低频天线，对钻井周围区域进行共中心点法(cmp法)测量。
步骤c.根据目标红壤关键带地下结构数据库，对目标红壤关键带雷达影像进行解译，获得雷达影像解译结果；其中，根据cmp探测推导相对介电常数，并计算不同层次结构的雷达波速；针对探地雷达数据解译：主要是使用各种算法对探地雷达数据进行处理，将红壤关键带地下结构的层位深度信息以图像的形式显示出来，如图7所示。
步骤d.基于雷达影像解译结果，更新目标红壤关键带地下结构数据库；其中，设定预测结果的空间分辨率为5m。
第五步、基于环境变量数据库构建三维预测模型，预测生成研究区不同层次结构厚度的空间分布图。
步骤e.基于目标红壤关键带地下结构数据库，结合环境变量数据库，构建目标红壤关键带地下结构三维预测模型，进而获得目标红壤关键带地下不同层次结构空间分布。
其中，根据验证样点集合的观测值与雷达影像解译结果对预测结果进行交叉验证，精度评价指标选用均方根误差。均质红土层、网纹红土层、风化层、总厚度的均方根误差分别为0.37m、0.67m、0.33m、0.72m，预测的空间分布图分别如图8、图9、图10、图11所示。
将所有层次厚度的预测结果进行空间叠加，生成新的图层，该图层表示了所有层次厚度的总和。将该图层与以depth_n为因变量生成的图层进行减操作，如图12所示，绝对值较高的地方表示了空间预测中不确定性较高。设计两条南北-东西方向的序列，如图13所示，可以得出不同序列地下结构层次的空间分布，如图14所示。最后，将生成的所有图层保存为栅格文件，用来表达关键带地下结构不同层次厚度的空间分布。
理解关键带各结构层次的厚度及其空间分布特征是气候变化、生态管护、水资源安全、自然灾害防治工作中日趋活跃的核心技术难题。所提出的基于探地雷达的红壤关键带地下结构空间预测方法不仅具有较好的迁移性，也能量化关键带地下结构的厚度与环境变量的关系。提出的预测方法具有“高内聚、低耦合”的技术特点，部分的技术细节能够较为方便地移植至其他相关的技术研发。所建立的一整套关键带层次结构的空间预测流程，有望为水资源安全、生态系统管护、气候变化等重大问题的解决提供了基础数据支持。
上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。
- 还没有人留言评论。精彩留言会获得点赞！