CN110780297B

CN110780297B - 一种采用星载合成孔径雷达监测土壤侵蚀的方法与系统

Info

Publication number: CN110780297B
Application number: CN201910942719.XA
Authority: CN
Inventors: 杨通; 张智杰; 郭旭东; 于潇; 韩圣其; 汪晓帆; 李潇; 刘智丽; 黄季夏; 南燕云
Original assignee: China Land Survey And Planning Institute
Current assignee: China Land Survey And Planning Institute
Priority date: 2019-09-30
Filing date: 2019-09-30
Publication date: 2023-05-23
Anticipated expiration: 2039-09-30
Also published as: CN110780297A

Abstract

本发明提供了一种采用星载合成孔径雷达监测土壤侵蚀的方法与系统。采用星载合成孔径雷达监测土壤侵蚀的方法包括：第一步骤，获取卫星平台的数据，并对获取的卫星平台的数据执行预处理；第二步骤，对预处理的卫星数据执行差分干涉处理，以生成差分干涉图；第三步骤，执行PS‑InSAR处理，生成PS点并迭代求解去除地形及大气相位；第四步骤，根据生成的差分干涉图获取累积沉降量；第五步骤，从累积沉降点中分离土壤侵蚀点，以土壤侵蚀点指示土壤侵蚀的空间分布及时间变化，并结合多源数据评估土壤侵蚀分布特征。本发明的采用星载合成孔径雷达监测土壤侵蚀的方法和系统的测量精度高，信噪比高，能够优化初始DEM，克服大气延迟，有利于土壤侵蚀的研究和管理。

Description

一种采用星载合成孔径雷达监测土壤侵蚀的方法与系统

技术领域

本发明涉及一种采用星载合成孔径雷达监测土壤侵蚀的方法与系统。

背景技术

土壤侵蚀是土壤及其母质在自然营力或人类活动的作用下被剥蚀、破坏、分离、搬运和沉积的过程。严重的土壤侵蚀会毁坏耕地资源、加剧泥沙下泄、淤积河流湖泊、影响防洪安全、诱发自然灾害、加重风沙扬尘、破坏自然和人文景观等，一方面迅速恶化生态环境，另一方面加剧当地人民群众的贫困程度。

现阶段的土壤侵蚀监测方法包括建立观测站、实地调查、立体摄影、高精度GPS、三维激光测量以及光学遥感方法等。遥感方法具有宏观性、综合性、可视性、实用性、经济性的特点，突破了传统的点线状土壤侵蚀微观机理研究范畴，将研究层面拓展到面上研究和动态研究中，因而对土壤侵蚀领域的研究发挥重要作用。

土壤侵蚀的主要特征之一是土层的剥蚀，中华人民共和国水利部发布的《土壤侵蚀分类分级标准SL190-2007》规定：水力侵蚀的土壤平均流失厚度介于5.9-11.1mm/a的，其侵蚀强度为极强烈；风力侵蚀的风蚀厚度介于10-25mm/a的，其侵蚀强度为中度。而合成孔径雷达的观测精度为厘米/毫米级，理论上满足中度以上土壤侵蚀监测的精度要求。

现阶段的土壤侵蚀遥感监测研究多基于光学遥感，监测的是正射水平面的二维特征，以像元流失作为土壤侵蚀的主要依据。相比之下，雷达遥感可以获取精准的高程变化信息，监测尺度由二维空间转换为三维空间，监测精度由米级平面精度转换为厘米/毫米级高程精度。

星载合成孔径雷达干涉测量InSAR(Synthetic Aperture RadarInterferometry)是雷达微波成像和电磁波干涉的集成技术，该技术通过合成孔径雷达在不同位置对观测区域的相位差等信息分析处理，进而获取观测区域的三维空间信息。InSAR技术以其覆盖范围大、不受云雾影响、厘米级监测精度的特点，成为监测地表形变最有前景的空间大地测量技术之一，并广泛应用于火山活动、地震、滑坡、不均匀沉降、地形测绘等研究领域。

但是，InSAR技术面临的困难之一，是植被和其他自然现象使地表反射率随时间发生变化，引起时间去相干，导致相位解缠难度增大而损失精度。故大多数传统的InSAR研究倾向于在工矿、道路、居民点等建成区开展。InSAR技术面临的困难之二，是雷达微波入射角不同，使斜距向地物频谱投影到数据频谱时发生偏移，导致空间去相干。InSAR技术面临的困难之三，是大气层使雷达微波发生速度变化和折射，导致大气延迟。

发明内容

本发明所要解决的技术问题是将InSAR技术引入土壤侵蚀监测中，根据识别目标的特点，优化某些技术环节，提供一种高测量精度、高信噪比的能够克服了大气延迟的采用星载合成孔径雷达监测土壤侵蚀的方法。

根据本发明，提供了一种采用星载合成孔径雷达监测土壤侵蚀的方法，包括：

第一步骤：获取卫星平台的数据，并对获取的卫星平台的数据执行预处理；

第二步骤：对预处理的卫星数据执行差分干涉处理，以生成差分干涉图；

第三步骤：执行PS-InSAR处理，生成PS点并迭代求解去除地形及大气相位；

第四步骤：根据生成的差分干涉图获取累积沉降量；

第五步骤：从累积沉降点中分离土壤侵蚀点，求解土壤侵蚀的空间分布及时间变化，并结合多源数据评估土壤侵蚀分布特征。

优选地，从累积沉降点中分离土壤侵蚀点包括：设置监测阈值以排除滑坡、泥石流等地表突变，执行相位标准差低通滤波以排除植被变化，执行历史形变量回溯以判定是否为土壤侵蚀、评估侵蚀营力。

优选地，根据生成的差分干涉图获取累积沉降点包括：从生成的差分干涉图中提取永久散射体，利用Delaunay三角网、星状网、树状网中的一种，根据提取的永久散射体来建立互相连结的永久散射体三角剖分网络，然后在去除残余相位因素之后获取累积沉降点。

优选地，采用振幅离差法和相干系数法的组合方法从生成的幅度和相干系数图中提取永久散射体。

优选地，去除残余相位因素分量包括：基于永久散射体三角剖分网络来计算永久散射体三角剖分网络上每条边的相邻永久散射体点的二次差分相位以得到二次差分相位模型，通过求解二次差分相位模型获得永久散射体点对的形变速率差和DEM修正差，通过相位解缠获取每个永久散射体点的绝对形变速率和DEM修正量，根据绝对形变速率和DEM修正量计算得到差分相位的形变相位和残余地形相位，执行滤波处理。

优选地，预处理包括：数据导入、精轨数据更新、主影像选择、DEM主辅影像配准、ESD配准、burst拼接、子条带拼接。

优选地，对预处理的卫星数据执行差分干涉处理包括如下步骤：首先对预处理的卫星数据执行ROI截取；随后针对截取的数据执行干涉相位计算处理，去除平地和地形相位，生成差分干涉图。

根据本发明，还提供了一种采用星载合成孔径雷达监测土壤侵蚀的系统，包括：

预处理单元，用于获取卫星平台的数据，并对获取的卫星平台的数据执行预处理；

差分干涉处理预处理单元，用于对预处理的卫星数据执行差分干涉处理，以生成差分干涉图；

PS-InSAR处理单元，用于生成PS点并迭代求解去除地形及大气相位；

累积沉降点获取单元，用于根据生成的差分干涉图获取累积沉降点。

处理单元，用于从累积沉降点中分离土壤侵蚀点，求解土壤侵蚀的空间分布及时间变化，并结合多源数据评估土壤侵蚀分布特征。

优选地，处理单元执行下述步骤以从累积沉降点中分离土壤侵蚀点：设置监测阈值以排除滑坡、泥石流等地表突变，执行相位标准差低通滤波以排除植被变化，执行历史形变量回溯以判定是否为土壤侵蚀、评估侵蚀营力。

优选地，根据生成的差分干涉图获取累积沉降量包括：从生成的差分干涉图中提取永久散射体，利用Delaunay三角网、星状网、树状网中的一种，根据提取的永久散射体来建立互相连结的永久散射体三角剖分网络，然后在去除残余相位因素之后获取累积沉降点。

根据本发明的采用星载合成孔径雷达监测土壤侵蚀的方法和系统的测量精度高，信噪比高，能够优化初始DEM，克服大气延迟，有利于土壤侵蚀的研究和管理。

附图说明

结合附图，并通过参考下面的详细描述，将会更容易地对本发明有更完整的理解并且更容易地理解其伴随的优点和特征，其中：

图1示意性地示出了根据本发明优选实施例的示例示出的研究区位置。

图2示意性地示出了根据本发明优选实施例的采用星载合成孔径雷达监测土壤侵蚀的方法的总体流程图。

图3示意性地示出了双基线图。

图4示意性地示出了根据本发明优选实施例的残余相位结果。

图5示意性地示出了根据本发明优选实施例的局部PS点分布图。

图6a至图6d示意性地示出了土壤侵蚀分布规律的示例。

图7示意性地示出了累计土壤侵蚀量的柱状图。

需要说明的是，附图用于说明本发明，而非限制本发明。注意，表示结构的附图可能并非按比例绘制。并且，附图中，相同或者类似的元件标有相同或者类似的标号。

具体实施方式

为了使本发明的内容更加清楚和易懂，下面结合具体实施例和附图对本发明的内容进行详细描述。

本发明拟利用卫星平台，使用星载永久散射体合成孔径雷达PS-InSAR(Persistent Scatterer Synthetic Aperture Radar Interferometry)技术，对内蒙古和林格尔县进行土壤侵蚀的面上观测，并对观测结果验证精度。同时，利用多源数据集成(光学遥感、降雨、坡度、土壤)，评估区域土壤侵蚀驱动因素及宏观分布规律。

星载永久散射体合成孔径雷达PS-InSAR克服了传统InSAR技术之不足，通过分析时间序列影像，提取时序内相位和幅度保持稳定的离散点，即PS(Persistent Scatterer，永久散射体)点，反演地表缓慢形变。，PS-InSAR聚焦于永久散射体(PS点)，削弱了农作物和森林的影响，将感兴趣区圈定在裸露地表，克服了植被区域去相干；PS-InSAR使用长时间序列点目标分析技术，通过形变速率估计、非线性形变和大气效应分离，克服了大气延迟。比之传统InSAR，PS-InSAR显著提高了测量精度和信噪比，高程向测量精度优于厘米级甚至达到毫米级。

<研究区示例>

研究区位于内蒙古自治区中南部(北纬39°58′－40°41′，东经111°26′－112°18′)，为首府呼和浩特市所辖旗县之一，蒙语“和林格尔”意为二十户人家。总面积3436平方公里，海拔1400～2028m，年平均降雨量为376mm，年均蒸发量为2073mm。和林格尔县属内蒙古高原和黄土高原的过渡地带，东南部山区、黄土丘陵区的面积占县域面积的77.7％；西北部属土默川平原的边缘，面积占和林格尔县总面积的22.3％。研究区超过一半区域属于黄土高原，黄土高原是世界上土壤侵蚀最严重区域，每年入黄泥沙量达16亿吨，且每年由于人类活动新增水土流失面积900-1100km2。自大规模农耕活动开始，人类活动就已经与水力、风力并列，成为土壤侵蚀的主要营力。根据《土壤侵蚀分类分级标准SL190-2007》分类体系，研究区既属于西北黄土高原水力侵蚀区，又属于内蒙古高原草原中度风蚀水蚀区，水力侵蚀与风力侵蚀并存。

在本申请说明书中作为示例示出的研究区位置见图1。

<数据源>

监测数据：2017至2018年的17期第一卫星平台(例如，哨兵一号Sentinel-1A)雷达数据，每期1景，用于监测土壤侵蚀。第一卫星平台的轨道类型为近圆形的太阳同步轨道，其轨道高度为693km，倾角为98.18°。卫星重访周期为12天，内有175个轨道。本研究所用的SAR数据是干涉宽模式(Interferometric Wide，IW)下的S1 TOPS-modes单视SLC(single lookcomplex)数据，SLC数据处理后的距离向分辨率为2.3m，方位向分辨率为13.9m，极化模式采用VV极化。

辅助数据：光学遥感为同步观测的第二卫星平台(哨兵二号Sentinel-2A)数据，用于目视验证监测结果及反演植被盖度；DEM(数字高程模型，Digital Elevation Model)数据用于控制地形影响及评估土壤侵蚀驱动因素；降雨量和土壤数据用于评估驱动因素。以雷达数据为基准数据集，辅助数据集尽量使用同步观测数据，若无同步观测数据则使用最近时相数据。数据使用情况详见表1。

表1数据源一览表

<方法实施例>

理论上，重复通过观测区域的卫星平台不遵循相同的路径，同时也会导致一些交叉道分离，若两次卫星平台通过时的目标高程特征发生偏移，则可计算出相关像素的干涉相位差。InSAR干涉测量的总干涉相位差是地形、大气、高程形变量及噪声的函数。为了监测土壤侵蚀引发的高程形变量，需要减去地形、大气等因素的贡献。PS-InSAR的优势在于可以优化初始DEM并估计大气影响，从而以极高的精度检测并绘制地表短期或长期的景观变化，具备监测土壤侵蚀的能力。

如图2所示，根据本发明优选实施例的采用星载合成孔径雷达监测土壤侵蚀的方法包括：

所述卫星平台的数据可以是作为主卫星平台的第一卫星平台的数据。

具体地，例如，预处理包括：数据导入、精轨数据更新、主影像选择、DEM主辅影像配准、ESD配准、burst拼接、子条带拼接等。

具体地，例如，对预处理的卫星数据执行差分干涉处理包括如下步骤：首先对预处理的卫星数据执行ROI(Region of interest，感兴趣区域)截取；随后针对截取的数据执行干涉相位计算处理，去除平地和地形相位，生成差分干涉图。

第四步骤：根据生成的差分干涉图获取累积沉降量；

在优选实施例中，从累积沉降点中分离土壤侵蚀点包括：设置监测阈值以排除滑坡、泥石流等地表突变，执行相位标准差低通滤波以排除植被变化，执行历史形变量回溯以判定是否为土壤侵蚀、评估侵蚀营力。

在优选实施例中，根据生成的差分干涉图获取累积沉降点包括：从生成的差分干涉图中提取永久散射体，利用Delaunay三角网、星状网、树状网中的一种或多种，并根据提取的永久散射体来建立互相连结的永久散射体三角剖分网络，然后在去除残余相位因素之后获取累积沉降点。

进一步优选地，采用振幅离差法和相干系数法的组合方法从生成的幅度和相干系数图中提取永久散射体。

进一步优选地，去除残余相位因素分量包括：基于永久散射体三角剖分网络来计算永久散射体三角剖分网络上每条边的相邻永久散射体点的二次差分相位以得到二次差分相位模型，通过求解二次差分相位模型获得永久散射体点对的形变速率差和DEM修正差，通过相位解缠获取每个永久散射体点的绝对形变速率和DEM修正量，根据绝对形变速率和DEM修正量计算得到差分相位的形变相位和残余地形相位，执行滤波处理。

下面对具体示例进行描述。

<雷达数据预处理>

第一卫星平台(主卫星平台)数据需要在极化模式筛选、SLC数据生成、精轨数据辅助矫正、DEM配准等的基础上，充分利用影像burst重叠区域信息，实现影像的精确配准。配准整幅影像后根据研究区选择burst及子条带，方便后续处理。主要预处理步骤见下。

数据导入：生成SLC(Single Looks Complex)单视复数据，同时利用各子条带及burst的重叠信息生成镶嵌的SLC图像和多视强度图(Multi Looks Intensity)。

轨道信息校正：利用ESA(European Space Agency)提供的精轨POD数据，修正轨道信息，提高查询列表精度、减小定时及轨道误差对干涉相位的影响。

主影像选取：从时序影像中根据时间和空间基线及多普勒质心频率，且成像时间处于非生长季的影像主影像，尽量使基线图呈星状发散，例如可以选择第8期影像作为主影像，主影像表征植被、农作物的生物量峰值期，首末影像表征植被、农作物的生物量低谷期，可以通过计算相位标准差，以低通滤波方式滤除茂密植被和农作物干扰，将感兴趣区定位于土壤裸露地区，排除植被变化干扰，具体参见图3。

外部DEM导入：将DEM与主影像像配准，得到每个像元的三维坐标信息，同时帮助后续去除地形相位引起的误差。

ESD配准：利用ESD方法生成主影像与副影像的几何关系(查询列表)，由于影像配准精度直接决定干涉测量精度，ESD配准方法可以利用burst重叠区域之间的相位差分信息，在强度互相关配准后进一步精确方位向偏移，以确保配准精度达到1/1000的要求。

Burst及子条带拼接：将分块存储的第一卫星平台数据根据研究区的位置提取出来，以便缩小处理数据量加快处理时间。

<差分干涉处理>

主影像分别与辅助影像生成17期差分干涉图，并结合卫星轨道参数、成像模型、参考DEM，去地形相位和平地相位。每一个干涉对都需要进行配准和重采样，形成干涉相位；利用轨道和外部的DEM去除平地效应和地形效应，形成时间序列的干涉相位。具体包括感兴趣区选择、干涉处理、地形相位去除、差分模型构建四个步骤。首先，将提取后的burst及子条带进行裁剪，进一步缩小数据量，提高后续处理速度。其次，将预处后的主副影像进行干涉处理，得到干涉图。然后，构建差分相位模型，求解每个像元的差分相位。最后，利用外部参考DEM数据和星历轨道参数，计算地形相位，并把地形相位和平地相位从干涉相位中去除，得到差分干涉相位图。主辅影像干涉相位可表示为公式(1)：

式(1)中

为地形相位；/>

为形变相位；/>

为大气延迟相位；/>

为平地相位；/>

为噪声相位。

<永久散射体时序处理>

提取永久散射体，并利用Delaunay三角网、星状网、树状网中的一种或多种，对识别的PS点建立互相连结的PS网络，然后去除残余相位及大气象相位的影响，最终得到累积沉降量。具体处理步骤见下。

PS候选点选取：选择对雷达波后向散射较强，并且在时序上较稳定的地物目标，如建筑物与岩石等。时序干涉处理是针对点的处理方式，所以选点是该方法的关键，这一过程可以将相干性差或者一些噪点去除，留下高相干且稳定的点进行分析，以提高监测精度。考虑到研究区面积较大、地形复杂、人工地表与自然地表同时存在，为保证PS候选点密度，本文采用两种选点方法，一是振幅离差法，二是相干系数法，以组合方法优化PS候选点的数量和质量，确保后续形变解算分析在参考点之间的传递。

构建三角剖分网：利用Delauany三角化方法构建三角剖分网络，以解决PS离散点的图像重构问题。同时，利用三角网络的冗余测量边，以区域网络平差方法控制形变测量误差，削弱大气相位影响，作为后续处理的基础。

计算二次差分相位：基于三角剖分网，计算网络上每条边的相邻PS点的二次差分相位。如前所述，大部分空间相关误差相位已经抵消，残余相位主要为高程和形变相位，通过求解二次差分相位模型，获得PS点对的形变速率差和DEM修正差。选取参考点，沿参考点积分可得每个PS点的绝对形变速率和DEM修正量，

相位解缠：可得差分相位的形变相位和残余地形相位。残余相位(图4)主要包括非线性形变和大气延迟相位，需通过滤波方法分离非线性形变和大气延迟。由于测量相位被限制在[-π,π]内，不能直接进行滤波，需要先对残余相位进行解缠，还原相位的模糊周期，才能对其作进一步的滤波处理。

大气滤波：大气延迟相位在时间域上表现为随机高频信号，在空间域上则表现为平滑低频信号；非线性形变在时间域上表现为低频信号，在频率域表现为低频信号；噪声相位则为时域高频，频域高频。因此，根据各相位成分的特点，可通过时空域高通、低通滤波分离各相位成分，从而得到形变相位，再由地理编码得出累计沉降量。

<辅助数据处理>

光学遥感数据：用于与雷达数据协同，目视验证精度，并计算NDVI(NormalizedDifference Vegetation Index,归一化差分植被指数)。对2017年春夏秋三个季节的第二卫星平台数据依次进行辐射定标、大气校正、数据镶嵌、数据裁剪、合成NDVI分布图。值得一提的是，第一卫星平台与第二卫星平台的系统校正参数相同，简化了几何校正和联合配准的处理步骤。

降雨量数据：该数据为研究区范围内19个气象观测站的同步观测数据，精确到小时。分别计算各观测站2017年3月至2018年3月的降雨量总和，然后根据观测站空间坐标，在Arcgis10.5平台下进行克吕金插值，生成雷达观测区间的降雨量分布图。

DEM数据：该数据来源于地理空间数据云下载的3景30米SRTM(AdvancedSpaceborne Thermal Emission and Reflection Radiometer Global DigitalElevation Model Version 2)数据。生成研究区坡度图。

土壤数据：在中科院南土所的1：100万土壤普查数据基础上，采用学界广泛接受的模型计算土壤可蚀性K因子，计算方法如公式(2)、(3)。

式中，K为土壤可蚀性因子，单位为t·h·(MJ·mm)-1；SAN砂粒含量(％)；SIL为粉粒含量(％)；CLA粘粒含量(％)；C为土壤有机碳含量，单位为g·kg-1。

<监测精度评估>

在县域内共计提取了多个PS点，点云密度满足监测精度需求。业界通常将年变化小于±1cm的区域视作稳定地表，而年沉降超过3cm的PS点极少，且可能属于滑坡范畴，不在本实施例的研究之列。基于此，提取全部年沉降介于1-3cm的PS点，以此沉降厚度对应《土壤侵蚀分类分级标准》的中度以上土壤侵蚀(风蚀、水蚀)。将满足上述要求的PS指示点(以下简称指示点)套合在2017年秋季的真彩色遥感影像上(第二卫星平台，例如Sentinel-2A)，采用窗口分析模式进行目视验证，如图5所示。

可知，指示点多分布于山坡、山谷区域，且该区域植被覆盖极少，具备土壤侵蚀发生的两大诱因：地形陡峭、植被稀疏。窗口南部山体的北坡发育有数条侵蚀沟，侵蚀沟从山脊延伸至山谷，符合典型的沟蚀特征。同时，注意到窗口的东部有大量植被覆盖，没有指示点分布；窗口西部山体较为平缓，指示点较少。总体而言，指示点分布符合肉眼认知及地理常识，侧面反映方法可行，结果稳定。

为评估土壤侵蚀监测精度，于2019年8月赴研究区开展了为期5天的野外精度验证，共涉及32处点云、103个PS点。实地验证结果表明，86.4％的PS点位于山坡及沟壑，符合典型的沟蚀特征，且大部分点位可见新近剥蚀出露的新鲜土层，PS点位与实际土壤侵蚀发生地较为吻合。监测误差的主要来源是误提了发生沉降的建筑物，因此，在应用PS-InSAR方法监测土壤侵蚀时，可以通过掩膜建筑用地提高监测精度。

<侵蚀分布规律>

PS-InSAR方法的机理是以面成点，一个PS点对应地面2.3*13.9m区域的平均情况。成像时以面成点，分析时则以点代面，以点云聚集区指示土壤侵蚀区。将年侵蚀量介于1-3cm的指示点云，套合在同步观测并计算的动态数据(NDVI分布图、降雨量插值图)及静态数据(坡度图、土壤可蚀性)上，可见土壤侵蚀区域的空间分布规律非常明显，见图6a至图6d。

总体而言，指示点云呈北东向45°分布，贯穿县域南北。此外，在县域西北、县域东南也呈局部聚集态势，这些区域均为潜在的中度以上土壤侵蚀区。

图6a为研究区2017年夏季的NDVI分布图，经专家判读可知，研究区的活跃耕地多分布于县域西北部，以及近北东向、近东西向两条河流附近；森林多分布于县域中东部。前已述及，茂密的森林和农作物是PS-InSAR方法的最大干扰源，故在计算相位标准差时以低通滤波方式去除，使雷达回波直接反映土壤信息。宏观上，指示点云多分布于黄土高原区域，该区域植被群落以草本和灌丛为主，且较为稀疏，光学遥感特征近似裸地，前人研究表明，植被盖度是土壤侵蚀的先导因子，对土壤侵蚀有明显的抑制作用。此外，在县域西北的集约化耕地区域，也见指示点云局部聚集区，结合光学遥感目视判读发现，指示点云位置多为建筑，说明此地建筑发生了不均匀沉降，推测为大规模农业用水引起地下水位下降所致。建筑沉降不属于本文研究范畴，但此现象应引起当地政府重视，采取必要的调查和预防措施。

图6b为县域内19个气象观测站的降雨量插值结果，可见县域降雨量分布极不平衡，整体呈东高西低分布，东部山区的年降雨量是西北部平原的3倍。一方面，水汽遇到山地会沿山坡抬升，气温逐步下降，水汽遇冷凝结形成降雨；另一方面，县域森林多分布于东部山区，森林的蒸腾作用会改善微气候，形成降雨。需要注意的是，县域西北极干旱地区，年均降雨量只有150mm左右，而此地区恰恰是该县的粮食主产区和集约化耕地分布区，且未见显著的河流、湖泊等地表水补给源，说明该区域的集约化农业大量依赖地下水，应注意地下水位下降引起的建筑沉降及生态恶化问题。结合指示点云分析，土壤侵蚀多发于县域西部干旱地区。说明相对于多雨气候，极端干旱气候减少了植被盖度，加剧了土壤侵蚀(Cuca B,Tzouvaras M,et al,2016)。

图6c为DEM数据生成的坡度图，少量平原分布于县域西北，其余则为大量山区丘陵地形。平原与山区分界线大致呈北东向45°延申，与指示点云总体分布特征平行。相应的，土壤侵蚀多分布于山区丘陵，平原地区相对较少。

图6d为土壤可蚀性分布图，由K因子计算得出，结合指示点云发现二者套合良好，规律明显。宏观上，指示点云聚集区多分布于土壤可蚀性较高区域，而在土壤可蚀性较低区域则分布较少且稀疏，遥感监测结果与模型计算结果相符。

县域降雨量分布极不平衡，整体呈东高西低分布，东部山区的年降雨量是西北部平原的3倍。一方面，水汽遇到山地会沿山坡抬升，气温逐步下降，水汽遇冷凝结形成降雨；另一方面，县域森林多分布于东部山区，森林的蒸腾作用会改善微气候，形成降雨。需要注意的是，县域西北极干旱地区，年均降雨量只有150mm左右，而此地区恰恰是该县的粮食主产区和集约化耕地分布区，且未见显著的河流、湖泊等地表水补给源，说明该区域的集约化农业大量依赖地下水，应注意地下水位下降引起的建筑沉降及生态恶化问题。结合指示点云分析，土壤侵蚀多发于县域西部干旱地区。说明相对于多雨气候，极端干旱气候减少了植被盖度，加剧了土壤侵蚀。

平原与山区分界线大致呈北东向45°延申，与指示点云总体分布特征平行。相应的，土壤侵蚀多分布于山区丘陵，平原地区相对较少。

宏观上，指示点云聚集区多分布于土壤可蚀性较高区域，而在土壤可蚀性较低区域则分布较少且稀疏，遥感监测结果与模型计算结果相符。

<侵蚀营力>

前人研究认为，土壤母质、土壤表面粗糙度、土壤表面坡度、植被覆盖度决定了土壤侵蚀速率，本研究亦佐证了这一结论。同时发现，在干旱半干旱地区，降雨量与土壤侵蚀是负相关的，难以判定土壤侵蚀的主要营力是风力抑或水力。PS-InSAR方法的另一个优势，在于时序数据支持对历史形变量的回溯。为厘清研究区土壤侵蚀的主要营力，对观测结果进行月度回溯，分别计算累计土壤侵蚀量与月度降雨量，见图7。其中，累计土壤侵蚀量为11161个土壤侵蚀点的平均侵蚀量，逐月累计。

回溯月度监测结果可知，研究区地表形变以缓慢变化的土壤侵蚀为主导，鲜见滑坡、泥石流等突发性形变。同时，研究区土壤侵蚀总量是逐月增加的，即便在没有降雨的12月、1月、2月，土壤侵蚀仍在发生，说明该区域土壤侵蚀的主要营力是风力，侵蚀类型以风蚀为主。此外，在降雨集中的6-9月，累计侵蚀量显著增加，变化梯度增大，说明降雨仍对研究区的土壤侵蚀具有驱动作用。据此推断，研究区土壤侵蚀的主要营力为风力，次要营力为水力，侵蚀类型为风蚀水蚀混合。11161个点年均侵蚀15mm，风蚀强度为中度。

<结论>

PS-InSAR方法理论上具备监测厘米/毫米级地表形变的能力，本发明通过设置监测阈值(高程减少1-3cm)、相位标准差低通滤波、历史形变量回溯等方式，排除滑坡、泥石流、植被变化等噪声影响，将监测目标锚定为土壤侵蚀。实证研究发现，该方法不但可以指示土壤侵蚀发生区域，还能得出土壤剥蚀的具体数值，反演侵蚀程度和侵蚀营力，在实现面上土壤侵蚀监测监管的同时，还可用于侵蚀驱动力建模。因此，该方法在土壤侵蚀研究领域具有广阔的应用前景。

研究区每年退耕还林一万亩，客观上抑制了土壤侵蚀的发生。在实地调查中发现，研究区实行针对性的退耕还林工程，大量的侵蚀沟附近均有人造林挡护，使土壤侵蚀至林线而止。然而，研究区大部属于黄土高原区，脆弱的自然本底决定了本区土壤侵蚀的多发、易发。此外，虽然陡坡开垦已被废止，但侵蚀沟旁的平坦耕地仍在种植，长此以往，势必造成“地越种越少、侵蚀愈演愈烈”的局面，直至此类耕地完全丧失农业价值和生态价值。因此，对于侵蚀沟旁的平坦耕地，应重视其生态敏感性，建议契合退耕还林、土地整治、第三次土地调查、生态保护红线等国家重大工程，变更其耕地性质，并以人工造林方式诱导其生态系统持续向好发展。

本发明将星载合成孔径雷达InSAR干涉测量技术应用于土壤侵蚀监测，探索雷达遥感监测宏观土壤侵蚀方法，证明了InSAR监测土壤侵蚀的可行性及应用前景。通过设置监测阈值、相位标准差低通滤波，从所有个PS点中分离出一部分沉降点，用以指示土壤侵蚀。经实地验证，提取的土壤侵蚀点与土壤侵蚀发生地的对应准确度为86.4％，年均侵蚀量为15mm。结合多源数据分析，研究区土壤侵蚀与植被盖度、降雨量、坡度、土壤成分密切相关。回溯时序观测结果发现，研究区土壤侵蚀的主要营力以风力为主，水力次之。该方法可为土壤侵蚀监测及定向提升侵蚀区管理提供助益。

<系统实施例>

可以理解的是，在本发明的另一些实施方式下，本发明还提供了一种采用星载合成孔径雷达监测土壤侵蚀的系统。

具体地，采用星载合成孔径雷达监测土壤侵蚀的系统可包括：

(1)预处理单元，用于获取卫星平台的数据，并对获取的卫星平台的数据执行预处理；

(2)差分干涉处理预处理单元，用于对预处理的卫星数据执行差分干涉处理，以生成差分干涉图；

(3)PS-InSAR处理单元，用于生成PS点并迭代求解去除地形及大气相位。

(4)累积沉降量获取单元，用于根据生成的差分干涉图获取累积沉降量。

(5)处理单元，用于从累积沉降点中分离土壤侵蚀点(设置监测阈值(排除滑坡、泥石流等地表突变)、相位标准差低通滤波(排除植被变化)，历史形变量回溯(判定是否为土壤侵蚀、评估侵蚀营力))，求解土壤侵蚀的空间分布及时间变化，并结合多源数据反演土壤侵蚀驱动力及分布规律。

在优选实施例中，根据生成的差分干涉图获取累积沉降量包括：从生成的差分干涉图中提取永久散射体，利用Delaunay三角网、星状网、树状网中的一种，并根据提取的永久散射体来建立互相连结的永久散射体三角剖分网络，然后在去除残余相位因素之后获取累积沉降点。

与其他领域的InSAR方法相比，所述方法的一些特点在于：1)使用了优化的PS点选点方法(幅度离差法+相干系数法)，增加了PS点的数量和密度，提高了监测精度；2)通过计算相位标准差，以低通滤波方式滤除茂密植被和农作物干扰，将感兴趣区定位于土壤裸露地区，排除植被变化干扰；3)通过设置监测阈值，排除了地形突变干扰，如滑坡、泥石流等；4)通过历史形变量回溯，判定是否为土壤侵蚀，及侵蚀营力。综上，从累积沉降点中，分离出土壤侵蚀点，达到监测土壤侵蚀的目的。

需要说明的是，除非特别指出，否则说明书中的术语“第一”、“第二”、“第三”等描述仅仅用于区分说明书中的各个组件、元素、步骤等，而不是用于表示各个组件、元素、步骤之间的逻辑关系或者顺序关系等。

可以理解的是，虽然本发明已以较佳实施例披露如上，然而上述实施例并非用以限定本发明。对于任何熟悉本领域的技术人员而言，在不脱离本发明技术方案范围情况下，都可利用上述揭示的技术内容对本发明技术方案作出许多可能的变动和修饰，或修改为等同变化的等效实施例。因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均仍属于本发明技术方案保护的范围内。

Claims

1.一种采用星载合成孔径雷达监测土壤侵蚀的方法，其特征在于包括：

第四步骤：根据生成的差分干涉图获取累积沉降点；

第五步骤：从累积沉降点中分离土壤侵蚀点，求解土壤侵蚀的空间分布及时间变化，并结合多源数据评估土壤侵蚀分布特征；

其中，从累积沉降点中分离土壤侵蚀点包括：设置监测阈值以排除包括滑坡、泥石流的地表突变，执行相位标准差低通滤波以排除植被变化，执行历史形变量回溯以判定是否为土壤侵蚀、评估侵蚀营力。

2.根据权利要求1所述的采用星载合成孔径雷达监测土壤侵蚀的方法，其特征在于，根据生成的差分干涉图获取累积沉降点包括：从生成的差分干涉图中提取永久散射体，利用Delaunay三角网、星状网、树状网中的一种，根据提取的永久散射体来建立互相连结的永久散射体三角剖分网络，然后在去除残余相位因素之后获取累积沉降量。

3.根据权利要求2所述的采用星载合成孔径雷达监测土壤侵蚀的方法，其特征在于，去除残余相位因素分量包括：基于永久散射体三角剖分网络来计算永久散射体三角剖分网络上每条边的相邻永久散射体点的二次差分相位以得到二次差分相位模型，通过求解二次差分相位模型获得永久散射体点对的形变速率差和DEM修正差，通过相位解缠获取每个永久散射体点的绝对形变速率和DEM修正量，根据绝对形变速率和DEM修正量计算得到差分相位的形变相位和残余地形相位，执行滤波处理。

4.根据权利要求1或2所述的采用星载合成孔径雷达监测土壤侵蚀的方法，其特征在于，预处理包括：数据导入、精轨数据更新、主影像选择、DEM主辅影像配准、ESD配准、burst拼接和子条带拼接。

5.根据权利要求1或2所述的采用星载合成孔径雷达监测土壤侵蚀的方法，其特征在于，对预处理的卫星数据执行差分干涉处理包括如下步骤：首先对预处理的卫星数据执行ROI截取；随后针对截取的数据执行干涉相位计算处理，去除平地和地形相位，生成差分干涉图。

6.一种采用星载合成孔径雷达监测土壤侵蚀的系统，其特征在于包括：

累积沉降点获取单元，用于根据生成的差分干涉图获取累积沉降点；

处理单元，用于从累积沉降点中分离土壤侵蚀点，求解土壤侵蚀的空间分布及时间变化，并结合多源数据评估土壤侵蚀分布特征；其中，处理单元执行下述步骤以从累积沉降点中分离土壤侵蚀点：设置监测阈值以排除包括滑坡、泥石流的地表突变，执行相位标准差低通滤波以排除植被变化，执行历史形变量回溯以判定是否为土壤侵蚀、评估侵蚀营力。

7.根据权利要求6所述的采用星载合成孔径雷达监测土壤侵蚀的系统，其特征在于，根据生成的差分干涉图获取累积沉降点包括：从生成的差分干涉图中提取永久散射体，利用Delaunay三角网、星状网、树状网中的一种，根据提取的永久散射体来建立互相连结的永久散射体三角剖分网络，然后在去除残余相位因素之后获取累积沉降量。

8.根据权利要求7所述的采用星载合成孔径雷达监测土壤侵蚀的系统，其特征在于，去除残余相位因素分量包括：基于永久散射体三角剖分网络来计算永久散射体三角剖分网络上每条边的相邻永久散射体点的二次差分相位以得到二次差分相位模型，通过求解二次差分相位模型获得永久散射体点对的形变速率差和DEM修正差，通过相位解缠获取每个永久散射体点的绝对形变速率和DEM修正量，根据绝对形变速率和DEM修正量计算得到差分相位的形变相位和残余地形相位，执行滤波处理。