CN113591572B - 基于多源数据和多时相数据的水土流失定量监测方法 - Google Patents

基于多源数据和多时相数据的水土流失定量监测方法 Download PDF

Info

Publication number
CN113591572B
CN113591572B CN202110722290.0A CN202110722290A CN113591572B CN 113591572 B CN113591572 B CN 113591572B CN 202110722290 A CN202110722290 A CN 202110722290A CN 113591572 B CN113591572 B CN 113591572B
Authority
CN
China
Prior art keywords
value
soil
water
factors
factor
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110722290.0A
Other languages
English (en)
Other versions
CN113591572A (zh
Inventor
沙晋明
徐赛萍
李小梅
吴龙滨
王金亮
包忠聪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fujian Normal University
Original Assignee
Fujian Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fujian Normal University filed Critical Fujian Normal University
Priority to CN202110722290.0A priority Critical patent/CN113591572B/zh
Publication of CN113591572A publication Critical patent/CN113591572A/zh
Application granted granted Critical
Publication of CN113591572B publication Critical patent/CN113591572B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0635—Risk analysis of enterprise or organisation activities
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A40/00—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production
    • Y02A40/10—Adaptation technologies in agriculture, forestry, livestock or agroalimentary production in agriculture
    • Y02A40/22—Improving land use; Improving water use or availability; Controlling erosion

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Strategic Management (AREA)
  • Economics (AREA)
  • Medical Informatics (AREA)
  • Development Economics (AREA)
  • Operations Research (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Game Theory and Decision Science (AREA)
  • Educational Administration (AREA)
  • General Business, Economics & Management (AREA)
  • Marketing (AREA)
  • Tourism & Hospitality (AREA)
  • Quality & Reliability (AREA)
  • Pit Excavations, Shoring, Fill Or Stabilisation Of Slopes (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种基于多源数据和多时相数据的水土流失定量监测方法，实施步骤如下：1)获取待监测区域的遥感影像，对原始遥感影像进行大气校正，根据COST模型的理论与方法，得到影像的地表反射率值；2)进行水土流失强度定性评价，其包括植被覆盖度、土地利用类型、沟壑密度和最优地形因子四个因子；3)利用遥感与GIS技术，结合多源数据，进行水土流失强度定量评价，根据USLE模型的计算公式A=R×K×LS×C×P，计算待监测区域的年均土壤流失量估值A并输出；本发明具有估算结果准确、通用性好的优点，能准确反应水土流失随时间变化的规律，实现区域水土流失的长期动态监测的优点。

Description

基于多源数据和多时相数据的水土流失定量监测方法
技术领域
本发明涉及水土流失风险评估领域，具体涉及一种基于USLE方程的水土流失遥感动态监测方法。
背景技术
水土流失是指由于水的侵蚀或者风力的作用使得土壤迁出土体，导致地力下降、严重的甚至完全失去地力。然而，要治理水土流失，我们必须明确水土流失的范围和程度。这样，才能有针对性地和开展有效地水土流失治理。
美国在二十世纪三十年代初发生严重的水土流失之后，成立了土壤保持局和水土保持国家实验室，开展了大量、细致的研究工作。经过30多年的努力，在1965年，得出了著名的通用水土流失方程(UniversalSoilLossEquation，USLE)。在此基础上，又经过近三十年的努力，获得了修订通用水土流失方程(ModifiedUniversalSoilLossEquation，MUSLE) 和修正通用水土流失方程(RevisedUniversalSoilLossEquation，RUSLE)。目前，国际上大多采用RUSLE来计算年均土壤流失量。其由六个因子决定：A＝R×K×LS×C×P，其中，A为年均土壤流失量估算值(tha-1yr-1)，R为降雨侵蚀力(MJmmha-1h-1yr-1)，K为土壤可蚀性因子(thahha-1MJ-1mm-1)，LS为坡长与坡度结合量(无量纲)，C为植被覆盖- 管理因子(无量纲)，P是水土保持措施因子(无量纲)。
根据以上公式，土壤侵蚀量由气候(降雨)、土壤、地形、植被和土地利用/覆盖等因子共同作用。其中，降雨、土壤和地形因子受自然条件影响，不同时期变化不大。而植被覆盖和土地利用等因子受人类活动影响，是易于变化但又是具有重要作用、影响极大的因子。
目前，在小范围或者局部田块，可以通过田间实测获得水土流失量。大范围的则通过遥感信息来获取，主要涉及两个方面：
1)直接利用遥感影像对水土流失或风险进行分级分类，包括目视解译和计算机自动解译。
目视解译一直是我国水土保持部门进行水土流失遥感调查采取的主要手段。我国于1985 年、1999年、2001年先后开展了全国第一、二、三次土壤侵蚀遥感调查均采用该方法。该方法的优点在于可以将人的经验和知识与遥感技术结合起来，充分利用专家的先验知识，避免了单纯的光谱分析可能带来的误差。然而解译没有明确的标准，解译过程中主观性极强，使得其结果难以在空间区域和时间序列上进行对比，而且其需要投入大量的人力、资金和时间，使得成本和时效不能兼顾。
计算机自动解译是为了避免采用目视解译技术所带来的大量人力、资金和时间的投入，从卫星影像中自动提取土壤侵蚀信息的另一种可行技术。但该方法基本上是依靠单纯的光谱信息进行分类，由于土壤侵蚀本身并不是以特定的土地覆盖等地表特征出现，而且指示土壤侵蚀的土壤属性光谱信息往往被植被覆盖、田间管理和耕种方式等等这样的土壤表层信息所掩盖，理论上只利用单独利用遥感信息是难以准确判断土壤侵蚀状况的。影像分类法在土壤侵蚀研究中的应用往往局限在某些特定的半干旱地区，而对于我国南方多云雨、地形复杂区域，限制了该方法的使用。
2)利用遥感数据提取影响水土流失的某些因子，然后利用一定的数据集成模式将因子综合得到水土流失或风险等级或分类图。
土壤侵蚀影响因子主要有降雨侵蚀力因子R，土壤可蚀性因子K，地形因子LS，植被覆盖-管理因子C和水土保持措施因子P。土壤侵蚀是这些因子共同作用的结果。在这些因子中，降水一般由气象数据获得，地形因子与土壤因子可以从遥感数据获取，但由于数据与技术的限制，目前还未见从遥感数据获取降雨侵蚀力因子的报道，地形数据更多是由等高线或DEM 数据产生，土壤因子更多从土壤图获取。所以，从遥感影像获得的主要是植被覆盖-管理因子和水土保持措施因子。也正是由于目前数据和技术的这种限制，利用遥感信息进行水土流失监测时，除了遥感数据，往往还要求有大量翔实的其它来源的数据辅助完成。但实际上，这些辅助数据的获取或更新存在一定难度，尤其在欠发达地区，这些数据在大面积上不易获得，要获取与遥感影像时相相匹配的数据更是难上加难。因此，当在一个区域或国家尺度时，数据的可获取性已成为水土流失动态监测的瓶颈之一。
由于数据的难获取性，现有方法要么只能在小范围通过田间实测获得土壤侵蚀量，若是在大范围，获得土壤侵蚀量不客观，和实际情况有较大差别，目前人们无法获得大范围的、有明确科学依据的水土流失信息，更无法进行区域水土流失的长期动态监测。
发明内容
本发明要解决的技术问题是提供一种基于多源数据和多时相数据的水土流失定量监测方法，其具有估算结果合理、通用性好、适合区域水土流失的长期动态监测、适用范围广的特点。
为解决上述技术问题，本发明采用的技术方案为：
一种基于多源数据和多时相数据的水土流失定量监测方法，其特征在于步骤如下：
1)获取待监测区域的原始遥感影像，对原始遥感影像进行大气校正，根据COST模型的理论与方法，从不同途径获取COST模型中所需要的参数值，按照COST模型计算流程对原始影像进行波段运算，得到影像的地表反射率值；
2) 进行水土流失强度定性评价，其包括植被覆盖度、土地利用类型、沟壑密度和最优地形因子四个因子；从步骤1）获得的原始遥感影像上提取出植被覆盖度；采用最大似然分类将监测区域的土地利用类型分为7类，得到土地利用类型图；对DEM进行处理，提取出沟壑分布情况，并与格网图叠加生成沟壑密度图；基于11个常用的地形因子，将其与坡长坡度因子做相关分析，选取相关系数最高的地形因子作为最优地形因子；采用主成分分析法计算这四个因子的权重值并利用GIS空间分析进行叠加，得到水土流失强度等级分布图；
3) 利用遥感与GIS技术，结合多源数据，进行水土流失强度定量评价，其坡长指数和坡度指数提取采用第一作者江忠善所提出的公式和经验值；植被覆盖与管理因子的提取是基于植被覆盖度，通过对第一作者蔡崇法所建立的模型进行修正后，按照修正后的模型公式计算得到；水土保持工程措施因子的获取是在土地利用类型的基础上，根据不同土地利用类型的工程措施情况进行赋值，根据USLE模型的计算公式A=R×K×LS×C×P，计算待监测区域的年均土壤流失量估值A并输出，其中，R为降雨侵蚀力因子，K为土壤可蚀性因子，LS为坡度坡长综合因子，C为植被覆盖和管理因子，P为水土保持措施因子。
上述步骤1)中根据COST模型的理论与方法，从不同途径获取COST模型中所需要的参数值，按照COST模型计算流程对原始遥感影像进行波段运算，得到原始遥感影像的地表反射率值，其具体步骤包括：
A1）将原始遥感影像DN值转换为光谱辐射值，即根据传感器的增益与偏移进行传感器辐射定标，根据获取辐射校正参数，其中/>为像元在传感器处的光谱辐射值，/>为以DN表示的经过量化定标的像元值，/>为8位DN值的理论最大值、取255，/>为根据/>拉伸的最大光谱辐射值，/>为根据8位DN值的理论最小值取0拉伸的最小光谱辐射值；
A2）计算原始遥感影像的大气程辐射值；
A3）进行大气校正，根据获取地球表面像元相对反射率，其中/>为/>波段中地表像元相对反射率；
所述步骤A2中计算原始遥感影像的大气程辐射值的详细步骤包括：
A21）根据得到传感器/>波段最小光谱辐射值，其中/>为传感器/>波段最小光谱辐射值，/>为/>波段最小DN值；
A22）根据得到假设黑体反射率为1%的/>波段黑体辐射值，其中/>指假设黑体反射率为1%的/>波段黑体辐射值，/>为/>波段的大气层顶平均太阳光谱辐照度，/>为太阳天顶角，/>为日地天文距离；
A23）根据得到大气程辐射值，其中/>为大气程辐射值，/>为传感器/>波段最小光谱辐射值，/>指假设黑体反射率为1%的/>波段黑体辐射值。
上述步骤2）中从所述原始遥感像上提取出植被覆盖度是根据像元二分模型得到，具体是指根据得到，其中FVC为植被覆盖度，NDVIveg为植被贡献的NDVI值，NDVIsoil为裸地贡献的NDVI值；
所述步骤2）中采用最大似然分类将监测区域的土地利用类型分为7类，得到土地利用类型图，具体是指将监测区域土地利用类型分为水体、耕地、低植被、高植被、裸地、建筑和道路7类；
所述步骤2）中对DEM进行处理，提取出沟壑分布情况，并与格网图叠加生成沟壑密度图，具体是指对原始DEM做填汁处理，生成无洼地DEM，通过对DEM每一个栅格点与它相邻的栅格点进行坡向分析，建立水流方向的数字矩阵，并进行栅格的追踪累积计算，建立汇流累积量的数字矩阵，设定阈值获得沟壑分布图，再利用GIS自动生成公里格网，把格网分布图与沟壑分布图叠加分析，统计出落在每个栅格内的沟的长度，除以单元格网面积即得到沟壑密度；
所述步骤2）中基于11个常用的地形因子，将其与坡长坡度因子做相关分析，选取相关系数最高的地形因子作为最优地形因子，具体步骤包括：
B1）根据DEM提取出坡度、坡向、坡度变率、坡向变率、剖面曲率、平面曲率、地形粗糙度、地形起伏度、高程变异系数、地表切割深度以及坡长这11个常用的地形因子；
B2）利用Arcmap软件随机生成样本点，在Excel软件中进行统计分析，得出其相关系数，作为评价的标准，以选出最适合用于评价水土流失的地形因子；
所述步骤2）中采用主成分分析法计算这四个因子的权重值并利用GIS空间分析进行叠加，得到水土流失强度等级分布图，其具体步骤包括：
C1）根据对各参评因子原始数据进行标准化处理，其中a为指标的标准化值，/>为指标值，/>分别为指标的最小值、最大值；
C2）根据对称矩阵计算相关系数矩阵；
C3）根据Jacobi可比法求出特征值与相应的单位特征向量；
C4）根据和/>计算主成分贡献率及累计贡献率，其中n为主成分总数，选取累计贡献率达85%~95%的特征值/>为所对应的第一、第二，…，第p个主成分代替原来的变量，其中p≤n，n为自然数；
C5）根据计算主成分载荷，其中lij为主成分载荷；
C6）根据得到各指标的权重值：其中lij为主成分载荷，/>为各指标的权重值。
上述步骤3）中坡长指数和坡度指数的提取采用第一作者江忠善所提出的公式和经验值，具体是指根据第一作者江忠善综合全国已有研究成果的基础上确立的全国统一的坡长指数取值：当θ≤5°时，m = 0.15；当 5°<θ≤12°时，m = 0.2；当12°<θ≤ 22°时，m =0.35；当 22°<θ<35°时，m = 0.45，同时，建立了不同坡度的坡长指数取值公式为，其中S0为使用百分比法表示的坡度，其数值为1.3~1.4；
所述步骤3）中植被覆盖与管理因子的提取是基于植被覆盖度，通过对第一作者蔡崇法建立的模型进行修正后，按照此模型公式计算得到，具体是指用回归分析方法建立了C与植被覆盖度fc以百分数表示之间的关系式，当fc大于9.632%时，将其带入公式求得的C值不大于1，具体根据
得到植被覆盖与管理因子，其中C为植被覆盖与管理因子，fc为植被覆盖度；
所述步骤3）中水土保持工程措施因子的获取，具体是指林地、草地、城镇、农村居民点以及未利用地般不采取任何保护措施，其P值为1；耕地确定其值的大小需要大量的观测小区实测资料，对采用某一专门措施前后的土壤流失量进行观测，分析观测数值得到不同水土保持措施下P值的大小。
本发明具有下述优点：
1、本发明基于COST模型的大气校正，在完成对遥感数据进行大气校正的同时也实现了地表反射率的反演，其优点是完全基于影像进行大气校正，不需要地表测量参数和大气参数，较为简单实用。
2、本发明通过宏观和微观角度，提取坡度、坡向等11个地形因子，通过分析这些因子与坡长坡度因子的相关度，选取与坡长坡度因子相关系数较高的地形因子用于水土流失强度评价，克服了一种地形因子很难准确具体地表达地形特征的问题。
3、本发明通过基于NDVI植被指数二分模型提取植被覆盖度信息，克服了由于大气影像地表湿度条件的改变以及地表粗糙度、土壤类型等的不同，NDVIsoil和NDVIveg会随着时间和空间而变化的问题，合理的选择了置信度5%，取相应的置信区间内NDVI的最小值和最大值分别作为NDVI soil和NDVIveg的值，提取的结果能较好的反应植被覆盖特征。
4、本发明提出了在无法获取年平均侵蚀模数的情况下，可以利用遥感影像数字处理的方法，结合地形数据，提取影响水土流失的主导因子：植被覆盖度、地形坡度、沟谷密度、土地利用类型，进行水土流失的研究，在技术上是完全可行的，且能较好的提取到影响水土流失的主导因子。
5、本发明在提取沟壑密度的实验过程中，根据前人的研究可知采用全流域汇流累积平均值作为阈值，提取的沟壑效果最佳，得到的沟壑密度能从宏观上真实反映地形破碎程度和地面被径流切割的程度，可作为水土流失等级划分的参考指标。
6、本发明通过反复试验发现当植被覆盖度fc大于9.632%时，利用第一作者蔡崇法建立的植被覆盖与管理因子与植被覆盖度之间的关系模型求得的C因子值才是0~1之间，其计算结果准确，较好的反映了地面植被覆盖对水土流失的综合作用。
7、本发明通过遥感和GIS技术，结合多源数据，提取出水土流失因子信息，从定性与定量两个角度分析水土流失状况，既能客观估算区域的水土流失量，具有估算结果准确、通用性好的优点，又能反应水土流失随时间变化的规律，实现区域水土流失的长期动态监测，可以广泛的应用于农业、环保、水利、国土等部门。
附图说明
图1为本发明实施例的基本流程示意图。
图2为本发明实施例得到的植被覆盖度图（左）和土地利用图（右）进行对照。
图3为本发明实施例得到的沟壑密度图（左）和坡长图（右）进行对照。
图4为本发明实施例的常用地形因子选取图。
图5为本发明实施例得到的坡度图（左）和坡向图（右）进行对照。
图6为本发明实施例选取的样本点分布图。
图7为本发明实施例得到的相关性分析图，其中：A)为坡度与坡长坡度因子关系图；B）为地形起伏度与坡长坡度因子关系图；C）为地形粗糙度与坡长坡度因子关系图；D）为地表切割深度与坡长坡度因子关系图；E）为坡度变率与坡长坡度因子关系图；F）为坡向与坡长坡度因子关系图；G）为坡长与坡长坡度因子关系图；H为高程变异系数与坡长坡度因子关系图；I）为坡向变率与坡长坡度因子关系图；J）平面曲率与坡长坡度因子关系图；K）为剖面曲率与坡长坡度因子关系图。
图8为本发明实施例得到的坡度等级图（左）和植被覆盖度等级图（右）进行对照。
图9为本发明实施例得到的沟壑密度等级图（左）和土地利用类型图（右）进行对照。
图10为本发明实施例得到的水土流失强度等级图。
图11为本发明实施例得到的降雨侵蚀力图（左）和土壤可蚀性图（右）进行对照。
图12为本发明实施例得到的水土保持措施因子（左）和水土流失量分级图（右）进行对照。
具体实施方式
如附图所示，本实施例基于多源数据和多时相数据的水土流失定量监测方法的实施步骤如下：
1)获取待监测区域的原始遥感影像，对原始遥感影像进行大气校正，根据COST模型的理论与方法，从不同途径获取COST模型中所需要的参数值，按照COST模型计算流程对原始影像进行波段运算，得到影像的地表反射率值；
2) 进行水土流失强度定性评价，其包括植被覆盖度、土地利用类型、沟壑密度和最优地形因子四个因子；从步骤1）获得的原始遥感影像上提取出植被覆盖度；采用最大似然分类将监测区域的土地利用类型分为7类，得到土地利用类型图；对DEM进行处理，提取出沟壑分布情况，并与格网图叠加生成沟壑密度图；基于11个常用的地形因子，将其与坡长坡度因子做相关分析，选取相关系数最高的地形因子作为最优地形因子；采用主成分分析法计算这四个因子的权重值并利用GIS空间分析进行叠加，得到水土流失强度等级分布图；
3) 利用遥感与GIS技术，结合多源数据，进行水土流失强度定量评价，其坡长指数和坡度指数提取采用第一作者江忠善所提出的公式和经验值；植被覆盖与管理因子的提取是基于植被覆盖度，通过对第一作者蔡崇法所建立的模型进行修正后，按照修正后的模型公式计算得到；水土保持工程措施因子的获取是在土地利用类型的基础上，根据不同土地利用类型的工程措施情况进行赋值，根据USLE模型的计算公式A=R×K×LS×C×P，计算待监测区域的年均土壤流失量估值A并输出，其中，R为降雨侵蚀力因子，K为土壤可蚀性因子，LS为坡度坡长综合因子，C为植被覆盖和管理因子，P为水土保持措施因子。
上述步骤1)中根据COST模型的理论与方法，从不同途径获取COST模型中所需要的参数值，按照COST模型计算流程对原始遥感影像进行波段运算，得到原始遥感影像的地表反射率值，其具体步骤包括：
A1）将原始遥感影像DN值转换为光谱辐射值，即根据传感器的增益与偏移进行传感器辐射定标，根据获取辐射校正参数，其中/>为像元在传感器处的光谱辐射值，/>为以DN表示的经过量化定标的像元值，/>为8位DN值的理论最大值、取255，/>为根据/>拉伸的最大光谱辐射值，/>为根据8位DN值的理论最小值取0拉伸的最小光谱辐射值；
A2）计算原始遥感影像的大气程辐射值；
A3）进行大气校正，根据获取地球表面像元相对反射率，其中/>为/>波段中地表像元相对反射率；
所述步骤A2中计算原始遥感影像的大气程辐射值的详细步骤包括：
A21）根据得到传感器/>波段最小光谱辐射值，其中/>为传感器/>波段最小光谱辐射值，/>为/>波段最小DN值；
A22）根据得到假设黑体反射率为1%的波段黑体辐射值，其中/>指假设黑体反射率为1%的/>波段黑体辐射值，/>为/>波段的大气层顶平均太阳光谱辐照度，/>为太阳天顶角，/>为日地天文距离；
A23）根据得到大气程辐射值，其中/>为大气程辐射值，为传感器/>波段最小光谱辐射值，/>指假设黑体反射率为1%的/>波段黑体辐射值。
上述步骤2）中从所述原始遥感像上提取出植被覆盖度是根据像元二分模型得到，具体是指根据得到，其中FVC为植被覆盖度，NDVIveg为植被贡献的NDVI值，NDVIsoil为裸地贡献的NDVI值；
所述步骤2）中采用最大似然分类将监测区域的土地利用类型分为7类，得到土地利用类型图，具体是指将监测区域土地利用类型分为水体、耕地、低植被、高植被、裸地、建筑和道路7类；
所述步骤2）中对DEM进行处理，提取出沟壑分布情况，并与格网图叠加生成沟壑密度图，具体是指对原始DEM做填汁处理，生成无洼地DEM，通过对DEM每一个栅格点与它相邻的栅格点进行坡向分析，建立水流方向的数字矩阵，并进行栅格的追踪累积计算，建立汇流累积量的数字矩阵，设定阈值获得沟壑分布图，再利用GIS自动生成公里格网，把格网分布图与沟壑分布图叠加分析，统计出落在每个栅格内的沟的长度，除以单元格网面积即得到沟壑密度；
所述步骤2）中基于11个常用的地形因子，将其与坡长坡度因子做相关分析，选取相关系数最高的地形因子作为最优地形因子，具体步骤包括：
B1）根据DEM提取出坡度、坡向、坡度变率、坡向变率、剖面曲率、平面曲率、地形粗糙度、地形起伏度、高程变异系数、地表切割深度以及坡长这11个常用的地形因子；
B2）利用Arcmap软件随机生成样本点，在Excel软件中进行统计分析，得出其相关系数，作为评价的标准，以选出最适合用于评价水土流失的地形因子；
所述步骤2）中采用主成分分析法计算这四个因子的权重值并利用GIS空间分析进行叠加，得到水土流失强度等级分布图，其具体步骤包括：
C1）根据对各参评因子原始数据进行标准化处理，其中a为指标的标准化值，/>为指标值，/>分别为指标的最小值、最大值；
C2）根据对称矩阵计算相关系数矩阵；
C3）根据Jacobi可比法求出特征值与相应的单位特征向量；
C4）根据和/>计算主成分贡献率及累计贡献率，其中n为主成分总数，选取累计贡献率达85%~95%的特征值/>为所对应的第一、第二，…，第p个主成分代替原来的变量，其中p≤n，n为自然数；
C5）根据计算主成分载荷，其中lij为主成分载荷；
C6）根据得到各指标的权重值：其中lij为主成分载荷，/>为各指标的权重值。
上述步骤3）中坡长指数和坡度指数的提取采用第一作者江忠善所提出的公式和经验值，具体是指根据第一作者江忠善综合全国已有研究成果的基础上确立的全国统一的坡长指数取值：当θ≤5°时，m = 0.15；当 5°<θ≤12°时，m = 0.2；当12°<θ≤ 22°时，m =0.35；当 22°<θ<35°时，m = 0.45，同时，建立了不同坡度的坡长指数取值公式为，其中S0为使用百分比法表示的坡度，其数值为1.3~1.4；
所述步骤3）中植被覆盖与管理因子的提取是基于植被覆盖度，通过对第一作者蔡崇法建立的模型进行修正后，按照此模型公式计算得到，具体是指用回归分析方法建立了C与植被覆盖度fc以百分数表示之间的关系式，当fc大于9.632%时，将其带入公式求得的C值不大于1，具体根据
得到植被覆盖与管理因子，其中C为植被覆盖与管理因子，fc为植被覆盖度；
所述步骤3）中水土保持工程措施因子的获取，具体是指林地、草地、城镇、农村居民点以及未利用地般不采取任何保护措施，其P值为1；耕地确定其值的大小需要大量的观测小区实测资料，对采用某一专门措施前后的土壤流失量进行观测，分析观测数值得到不同水土保持措施下P值的大小。
本发明以平潭为具体实施例来证明本发明的方法是可行的：
一、水土流失强度定性评价：
本发明实施例中使用的多源数据为由水土流失项目提供的2003年5米分辨率DEM影像、土壤类型图（含K值）、降雨侵蚀力因子图和平潭乡镇行政边界矢量图以及由美国地质调查局网站（http://glovis.usgs.gov）提供的2009年Landsat5 TM原始遥感影像，影像的获取时间是2009年11月8日，轨道号是118/42，空间分辨率为30m（其中，热红外波段的空间分辨率为120m），太阳高度角约为63.39。。
1指标信息提取
水土流失主要受降雨、土壤母质、地形、植被覆盖、水土保持措施等因子的影响，是自然因素和人为因素综合作用的结果。据研究，当地面坡度在0。~40。情况下，坡度与土壤冲刷量呈正比关系；但水土流失量的大小，还与地表有无植被覆盖有密切的关系，大量实验结果表明，地形因素对水土流失的影响作用是受植被制约的，只有在植被被破坏后，地形才是造成侵蚀的重要因素。由于降雨侵蚀力和土壤可蚀性在较短的时间内不会发生大的变化，因此对水土流失强度等级评价可以基于植被覆盖度、地形和土地利用类型等因子进行。由于一种地形因子很难准确具体地表达地形特征，因此，有必要探讨最适宜地形因子的选取。本发明根据监测区域的实际地理情况，从宏观和微观角度，提取坡度、坡向等11个地形因子，通过分析这些因子与坡长坡度因子的相关度，选取与坡长坡度因子相关系数较高的地形因子用于水土流失强度评价。
1.1植被覆盖度
植被覆盖度(FVC,Fractional Vegetation Cover)是植被在地面的垂直投影面积占总面积的百分比，是侵蚀动力的抑制因子。然而，植被盖度与水土流失强度的关系并非简单的线性关系。换言之，植被覆盖度越高，则地表的抗侵蚀性越强。植被指数和植被覆盖度之间具有较高的线性正相关特性，像元二分模型是近年来发展起来的常用于由植被指数来确定植被覆盖度的一种有效的方法。通过计算得到的植被覆盖图如图2左边图所示，影像值为0到1之间。其中，颜色越浅的区域植被覆盖度越大，反之，颜色越深的区域植被覆盖度越小。
1.2土地利用类型
人类不合理的土地利用是水土流失发生的主要原因之一，因此，土地利用作为人类利用土地各种活动的综合反映和水土流失有着密切的联系。本发明采用监督分类中的最大似然分类方法，将监测区域土地利用类型分为水体、耕地、低植被、高植被、裸地、建筑和道路7类，如图2右边图所示。利用ENVI的统计功能可以得到各个地类占监测区域的面积比例。由表1可知，监测区域中大部分是低植被覆盖区，其次是耕地。
表1 土地利用类型统计
地类 像元数 面积比例
水体 24581 8.21%
耕地 55994 18.69%
低植被 130538 43.57%
高植被 38167 12.74%
裸地 16946 5.66%
建筑 11881 3.97%
道路 21478 7.17%
1.3沟壑密度
沟壑密度也称沟谷密度或沟道密度，指单位面积内沟壑的总长度。沟壑的形成过程是流域地貌演化的过程，也是水土流失的过程，因此，沟壑密度能从宏观上真实反映地形破碎程度和地面被径流切割的程度，可作为水土流失等级划分的参考指标。沟壑密度的数学表达式如下，其值越大，地表与降雨、径流接触发生侵蚀的面积就越大，降水径流的冲刷力和侵蚀力就越大，地表物质稳定性降低，越容易引发水土流失。
式中，Ds指沟壑密度（单位：公里/公里2），指样区内的沟壑总长度（单位：公里）；A指特定样区的面积（单位：公里2）。
由上述公式可知，为了得到沟壑密度，需要先提取出沟壑信息。具体提取方法：首先对原始DEM做填汁处理(Fill sink),生成无洼地DEM，通过对DEM每一个栅格点与它相邻的栅格点进行坡向分析,建立水流方向的数字矩阵(Flow direetion)，并进行栅格的追踪累积计算,建立汇流累积量的数字矩阵(Flow accumulation)，设定阈值获得沟壑分布图。在沟壑提取过程中，DEM网格尺寸和汇流累积阈值这两个因素起到重要作用。本文DEM分辨率为5m，已经是高分辨率影像，对于汇流累积阈值的取值，采用全流域汇流累积平均值作为阈值，提取的沟壑效果最佳。得到沟壑分布图之后再利用 GIS 自动生成监测区域公里格网，把格网分布图与沟壑分布图叠加分析，统计出落在每个栅格内的沟的长度，除以单元格网面积即得到沟壑密度，经整饰输出沟壑密度图，如图3的左边图。
1.4最优地形因子的选取
GIS数字地形分析常用因子分析法，不同的地形因子从不同侧面反映了地形特征。从其所描述的空间区域范围，常用的地形因子可以划分为微观地形因子和宏观地形因子两种基本类型[98]，常用的地形因子分类体系见图4。
常用地形因子的提取
（1）坡长
坡长通常是指在地面上一点沿水流方向到其流向起点间的最大地面距离在水平面上的投影长度。其数学表达为：
式中L指坡长，m指地表面沿流向的水流长度，θ指水流地区的地面坡度值。坡长图见图3的右边图。
（2）坡度和坡向
地表面任一点的坡度是指过该点的切平面与水平地面的夹角。坡度是构成水力侵蚀的动力条件，一般情况下，随着坡度的增加，水土流失强度会加剧，因此，坡度的大小是决定水土流失强弱的重要因素。坡度见图5的左边图被广泛应用于径流流速、植被、降雨量、地貌和土地适宜性评价等。坡向是指地表面上一点的切平面的法线矢量在水平面的投影与过该点的正北方向的夹角。对于地面任何一点来说，坡向见图5的右边图表征了该点高程值改变量的最大变化方向。不同坡向植被覆盖度不同，植物类型不同，对水土流失影响的程度不同。坡度和坡向是相互联系的两个地形参数，均可利用ArcGIS软件表面分析功能提取。对于坡向数据，坡向值有以下规定: 正北方向为0°，正东方向为90°，正南方向为180°，正西方向为270°。
（3）坡度变率和坡向变率
坡度变率，是地面坡度在微分空间的变化率，是依据坡度的求算原理，在所提取的坡度值的基础上对地面每一点再求算一次坡度。即坡度之坡度（Slope of Slope, 简称SOS）。坡度是地面高程的变化率的求解，因此，坡度变率表征了地表面高程相对于水平面变化的二阶导数。
坡向变率，是指在地表的坡向提取基础之上，进行对坡向变化率值的二次提取，亦即坡向之坡度（Slope of Aspect, SOA）。它可以很好的反映等高线弯曲程度。
（4）地面曲率
地面曲率是对地形表面一点扭曲变化程度的定量化度量因子，其在垂直和水平两个方向上的分量分别称为剖面曲率（Profile Curve）和平面曲率（Plan Curve）。剖面曲率是对地面坡度的沿最大坡降方向地面高程变化率的度量。平面曲率描述的是地表曲面沿水平方向的弯曲、变化情况。提取方法：利用ArcGIS软件中Spatial Analyst Tools工具集下的Surface子工具集里的Curvature工具实现。
（5）地形粗糙度和地形起伏度
地形粗糙度是指在一个特定的区域内，地球表面积与其投影面积之比。它也是反映地表形态的一个宏观指标。计算地形粗糙度的数学表达式为：
式中，angle表示地表面与水平面的坡度；R是地表粗糙度。由于cos（angle）用弧度制作为单位，因此，需要将所提取的坡度图中的角度转化为弧度。地形粗糙度是能反映地形的起伏变化和侵蚀程度的宏观地形因子。在研究水土流失监测具有重要意义。
地形起伏度是指在一个特定的区域内，所有栅格中最大高程与最小高程的差值。它是描述一个区域地形的一个宏观性的指标。提取方法：首先利用Spatial Analysis下使用栅格邻域计算工具求出一定范围内海拔高度的最大值和最小值，再对其求差值即可。计算公式为：
QFD=Hmax-Hmin
式中，QFD是分析区域内的地形起伏度，Hmax和Hmin分别是分析窗口内的最大高程值和最小高程值。在水土流失研究中，地形起伏度指标能够反映水土流失类型区的水土流失特征，是比较适合估计区域水土流失评价的地形指标。
（6）高程变异系数和地表切割深度
高程变异是反映分析区域内地表单元格网各顶点高程变化的指标，它以格网单元顶点的标准差s与平均高程z的比值来表示：/>
式中，V为高程变异系数，n是单元格网顶点数目，s为标准差。
地表切割深度是指地面某点的邻域围的平均高程与该邻域范围内的最小高程的差值。地表切割深度直观的反映了地表被侵蚀切割的情况并对这一地学现象进行了量化，是研究水土流失及地表侵蚀发育状况时的重要参考指标。其计算公式为 ：
式中，Di指地面每一点的地表切割深度，指一个固定分析窗口内的平均高程，指一个固定分析窗口内的最低高程。
地形因子与坡长坡度因子的相关性分析
USLE模型中的坡长坡度因子LS的提取方法与因子图祥。为了选取具有代表性的地形因子，本文将以上11个地形因子与LS因子做相关分析。利用Arcmap软件随机生成监测区域内15000个样本点（样本点分布见图6），样本点均匀布满监测区域范围，可靠性较高。分别导出LS与这11个地形因子的样本值，在Excel软件中进行统计分析。从LS与这11个地形因子的散点图（如图7）可以看出地形因子随着LS所变化的趋势，进而得出其相关系数,作为评价的标准,以选出最适合用于评价水土流失的地形因子。
图7汇总了LS与11个地形因子的相关性，从图上可以很直观地看出，相关系数最高者为坡度，其次是地表起伏度。LS与地表切割深度、地形粗糙度的相关系数也高于0.9。相关系数越高，表明该地形因子和水土流失地形分量的越接近，即其对水土流失的影响作用越明显。由此，可以初步选定坡度作为评价水土流失的基本地形因子，地表起伏度、地表切割深度、地形粗糙度为备选因子。
2主成分分析法
主成分分析是把多个变量划为少数几个综合指标的一种统计分析方法，从数学角度来看，这是一种降维处理技术。利用主成分分析法计算因子权重的步骤归纳如下：
（1）各参评因子原始数据的标准化处理：为实现各种评价参数之间的可比较性，需要对评价参数指标值进行标准化处理。
式中，a为指标的标准化值；为指标值；/>分别为指标的最小值、最大值。
（2）计算相关系数矩阵
本文共6个评价因子，因此其相关系数矩阵为6阶实对称矩阵R，公式如下：
/>
式中， (i，j=1，2，…，6)是第i个参评因子和第j个参评因子之间的相关系数；n为监测区域像元点总数，/>分别为第i个因子影像和第j个因子影像的第k个像元点的DN值。因为R是实对称矩阵，即/>，故只须要计算上三角元素或着下三角元素即可。
（3）计算特征值与相应的单位特征向量
解特征方程，通常用雅可比法（Jacobi）求出特征值/>，然后分别求出对应于/>的单位特征向量/>表示向量/>的第j个分量。
（4）计算主成分贡献率及累计贡献率
第i个主成分的贡献率和主成分累计贡献率/>的计算公式为：
式中，n为主成分总数。一般选取累计贡献率达85%~95%的特征值所对应的第一、第二，…，第p（p≤n）个主成分代替原来的变量。这前p个主成分称为共因子。
（5）计算主成分载荷
主成分载荷的计算公式为
（6）计算各指标的权重值
第j个变量在全部p个公共因子上的载荷的平方和叫做变量的公共性。/>的大小反映了变量j在公共性部分中的作用或重要性程度。因此，可以把每个变量对应的/>作为该变量的权重。计算公式为：
本发明根据主成分分析法计算得到的各个评价参数的权重值见表2。
表2 各个评价参数的权重值
权重值 植被覆盖度 沟壑密度 土地利用类型 坡度
hj 0.576 0.121 0.218 0.162
3水土流失空间分布分析
参考 1985 年全国水土流失遥感调查工作技术细则，根据《水土保持技术规范》的有关要求，确定坡度、植被覆盖度、沟壑密度和土地利用类型的定级标准划分等级（表3）并对这四个因子影像进行重分类，结果见图8和图9。
表3 各指标定级标准
微度 轻度 中度 强烈 极强烈 剧烈
坡度（°） <5 5-8 8-15 15-25 25-35 >35
植被覆盖度（%） >90 70-90 50-70 30-50 10-30 <10
沟壑密度（km/km2） 0-1 1-2 2-3 3-5 5-7 >7
土地利用类型 高植被、水体 中低植被 草地 建筑、道路 裸地 耕地
由植被覆盖度等级图可以看出，中低覆盖度主要分布于人类活动较为频繁的地区，如建筑、道路等，地形较为破碎，沟谷密集。中高覆盖多见于山体中上部，植被茂密。
由图9可知，平潭的水土流失主要是微度等级，其次是剧烈等级。水土流失出现恶化情况，主要是由于人类活动的不断增加，如过度放牧、乱砍滥伐等。监测区域水土流失分布与植被发育程度相关性很大。监测区域水土流失较为严重的区域主要分布于植被覆盖不良的荒地和坡耕地上，其中坡耕地水土流失最为严重，大多为强烈流失。从地貌上看，水土流失主要发生在沟谷和山坡中下部，这些地方人类活动频繁，植被破坏严重，大量开垦坡耕地，导致了水土流失加剧。该地区多为农耕地，土地的保护对农业生产的发展关系极大，因此应大力加强对人为破坏的管制，减少放牧、农耕对地表植被的破坏，并加强对区内水土流失的治理。
二、水土流失定量评价：
雨滴击溅作用和因降雨产生的径流，是土壤水蚀的动力因子。一般认为，在同等情况下，降水强度越大，水土流失强度越强烈。降雨侵蚀力因子，指降雨动能与降雨时段最大雨强的乘积，是评价降雨对土壤剥离、搬运侵蚀的动力指标，反映为当其它因子相同时，不同的气候和降雨对水土流失的作用，与降雨量、降雨历时、降雨强度、降雨动能有关。本研究中的降雨侵蚀力因子（R）图如图11左边图。
2.1土壤可蚀性因子
土壤可蚀性因子K反映土壤自身特性，是对降雨的冲刷和溅蚀作用抵抗能力的表述，定义为单位面积上每降雨侵蚀系数单位的土壤流失率。土壤可蚀性与土壤的理化性质密切相关，而土壤的理化性质受到成土母质和成土过程的影响，因此对于 K 值的描述需要综合多种特性，为此土壤学家做了大量的研究，提出多种计算方法，主要有公式法、查图表法、查诺谟法、模拟降雨实测法和自然降雨实测法。本文的土壤可蚀性因子见图11右边图。
2.2坡度坡长因子
地形是通过影响地表径流特征作用于水土流失，其中最为重要的因子是坡度和坡长因子。在USLE模型中，坡长坡度因子是在相同条件下，每单位面积坡面流失与标准小区(坡长22.1m，坡度9％)流失之比值。坡长坡度LS 因子的计算公式为：
式中，λ为坡长，单位为米；θ为坡度，单位为度；m为坡长指数；n为坡度指数。第一作者江忠善综合全国已有研究成果的基础上确立的全国统一的坡长指数取值：当θ≤5°时, m= 0.15;当 5°<θ≤12°时, m = 0.2;当12°<θ≤ 22°时, m = 0.35;当 22°<θ<35°时,m =0.45。同时,建立了不同坡度的坡长指数取值公式为，S0为使用百分比法表示的坡度，并推荐全国统一采用的坡度指数值为可取1.3~1.4，本文中n值取值为1.35。/>
2.3植被覆盖与管理因子
第一作者蔡崇法通过径流小区的多场人工降雨和天然降雨资料与地表覆盖度之间的相关关系，用回归分析方法建立了C与植被覆盖度fc (以百分数表示)之间的关系式。本发明通过反复试验可知当fc大于9.632%时，将其带入公式求得的C值才不大于1。因此本文对公式（1）进行必要的修改，修订模型见公式（2）。
2.4水土保持工程措施因子
根据土地利用现状及相关文献进行P值设定，见表4。生成P值栅格数据。
表4 平潭县典型土地利用类型P值
土地类型 水体 耕地 中低植被 高植被 裸地 道路 建筑
P值 0 0.35 0.4 0.2 1 0 0
3水土流失量分析
利用Arcmap软件的栅格计算器功能，将上述获取的各因子图层进行连乘。由于USLE 模型的单位是英制，所以应当将其计算结果乘以224.2 转换为公制单位 t·hm-2·a-1，得到各像元的土壤流失量。根据水土流失的试验结果，参考中国水利部 2007 年制定的《水土流失分类分级标准》，将水土流失强度等级划分为微度、轻度、中度、强度、极强烈和剧烈6个等级（表5），等级分布图见图12的右边图。
表5 水土流失强度分级标准
级别 平均侵蚀模数（t·km2）
微度 <1000
轻度 1000~2500
中度 2500~5000
强度 5000~8000
极强烈 8000~15000
剧烈 ＞15000
从水土流失量分级分布图（见图12）中可以看出，水土流失主要分布在海坛岛的东南部，即人为活动活跃区。植被覆盖度高的区域由于植被对侵蚀动力的抑制作用，不易发生水土流失。轻微流失区在监测区域内分布广泛，面积最大。该等级区的植被覆盖度较高，林地立地条件较好，土壤较为肥沃，坡度虽然也很大，但是岩石的抗侵蚀能力较强。在建筑区等人为活动剧烈区，植被遭到严重破坏，沟谷较发育，坡耕地上由于降雨侵蚀，常有细沟发育，沟谷较为密集，水土流失较为严重。
以上所述仅为本发明的优选实施方式，本发明的保护范围并不仅限于上述实施方式，凡是属于本发明原理的技术方案均属于本发明的保护范围。对于本领域的技术人员而言，在不脱离本发明的原理的前提下进行的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (1)

1.一种基于多源数据和多时相数据的水土流失定量监测方法，其特征在于步骤如下：
1)获取待监测区域的原始遥感影像，对原始遥感影像进行大气校正，根据COST模型的理论与方法，从不同途径获取COST模型中所需要的参数值，按照COST模型计算流程对原始影像进行波段运算，得到影像的地表反射率值；
2) 进行水土流失强度定性评价，其包括植被覆盖度、土地利用类型、沟壑密度和最优地形因子四个因子；从步骤1）获得的原始遥感影像上提取出植被覆盖度；采用最大似然分类将监测区域的土地利用类型分为7类，得到土地利用类型图；对DEM进行处理，提取出沟壑分布情况，并与格网图叠加生成沟壑密度图；基于11个常用的地形因子，将其与坡长坡度因子做相关分析，选取相关系数最高的地形因子作为最优地形因子；采用主成分分析法计算这四个因子的权重值并利用GIS空间分析进行叠加，得到水土流失强度等级分布图；
3) 利用遥感与GIS技术，结合多源数据，进行水土流失强度定量评价，其坡长指数和坡度指数提取采用以下所提出的公式和经验值，具体是指根据综合全国已有研究成果的基础上确立的全国统一的坡长指数取值：当θ≤5°时，m = 0.15；当 5°<θ≤12°时，m = 0.2；当12°<θ≤ 22°时，m = 0.35；当 22°<θ<35°时，m = 0.45，同时，建立了不同坡度的坡长指数取值公式为 ，其中S0为使用百分比法表示的坡度，其数值为1.3~1.4；植被覆盖与管理因子的提取是基于植被覆盖度，通过对建立的模型进行修正后，按照修正后的模型公式计算得到；水土保持工程措施因子的获取是在土地利用类型的基础上，根据不同土地利用类型的工程措施情况进行赋值，根据USLE模型的计算公式A=R×K×LS×C×P，计算待监测区域的年均土壤流失量估值A并输出，其中，R为降雨侵蚀力因子，K为土壤可蚀性因子，LS为坡度坡长综合因子，C为植被覆盖和管理因子，P为水土保持措施因子；
所述步骤1)中根据COST模型的理论与方法，从不同途径获取COST模型中所需要的参数值，按照COST模型计算流程对原始遥感影像进行波段运算，得到原始遥感影像的地表反射率值，其具体步骤包括：
A1）将原始遥感影像DN值转换为光谱辐射值，即根据传感器的增益与偏移进行传感器辐射定标，根据 获取辐射校正参数，其中/>为像元在传感器处的光谱辐射值， />为以DN表示的经过量化定标的像元值，/>为8位DN值的理论最大值、取255，/>为根据/>拉伸的最大光谱辐射值，/>为根据8位DN值的理论最小值取0拉伸的最小光谱辐射值；
A2）计算原始遥感影像的大气程辐射值；
A3）进行大气校正，根据获取地球表面像元相对反射率，其中/>为/>波段中地表像元相对反射率；
所述步骤A2中计算原始遥感影像的大气程辐射值的详细步骤包括：
A21）根据得到传感器/>波段最小光谱辐射值，其中/>为传感器/>波段最小光谱辐射值，/>为/>波段最小DN值；
A22）根据得到假设黑体反射率为1%的/>波段黑体辐射值，其中/>指假设黑体反射率为1%的/>波段黑体辐射值，/>为/>波段的大气层顶平均太阳光谱辐照度，/>为太阳天顶角，/>为日地天文距离；
A23）根据得到大气程辐射值，其中/>为大气程辐射值，为传感器/>波段最小光谱辐射值，/>指假设黑体反射率为1%的/>波段黑体辐射值；
所述步骤2）中从所述原始遥感像上提取出植被覆盖度是根据像元二分模型得到，具体是指根据得到，其中FVC为植被覆盖度，NDVIveg为植被贡献的NDVI值，NDVIsoil为裸地贡献的NDVI值，其中NDVI为归一化植被指数；
所述步骤2）中采用最大似然分类将监测区域的土地利用类型分为7类，得到土地利用类型图，具体是指将监测区域土地利用类型分为水体、耕地、低植被、高植被、裸地、建筑和道路7类；
所述步骤2）中对DEM进行处理，提取出沟壑分布情况，并与格网图叠加生成沟壑密度图，具体是指对原始DEM做填汁处理，生成无洼地DEM，通过对DEM每一个栅格点与它相邻的栅格点进行坡向分析，建立水流方向的数字矩阵，并进行栅格的追踪累积计算，建立汇流累积量的数字矩阵，设定阈值获得沟壑分布图，再利用GIS自动生成公里格网，把格网分布图与沟壑分布图叠加分析，统计出落在每个栅格内的沟的长度，除以单元格网面积即得到沟壑密度；
所述步骤2）中基于11个常用的地形因子，将其与坡长坡度因子做相关分析，选取相关系数最高的地形因子作为最优地形因子，具体步骤包括：
B1）根据DEM提取出坡度、坡向、坡度变率、坡向变率、剖面曲率、平面曲率、地形粗糙度、地形起伏度、高程变异系数、地表切割深度以及坡长这11个常用的地形因子；
B2）利用Arcmap软件随机生成样本点，在Excel软件中进行统计分析，得出其相关系数，作为评价的标准，以选出最适合用于评价水土流失的地形因子；
所述步骤2）中采用主成分分析法计算这四个因子的权重值并利用GIS空间分析进行叠加，得到水土流失强度等级分布图，其具体步骤包括：
C1）根据对各参评因子原始数据进行标准化处理，其中a为指标的标准化值，/>为指标值，/>分别为指标的最小值、最大值；
C2）根据对称矩阵计算相关系数矩阵；
C3）根据Jacobi可比法求出特征值与相应的单位特征向量；
C4）根据和/>计算主成分贡献率及累计贡献率，其中n为主成分总数，选取累计贡献率达85%~95%的特征值/>为所对应的第1、第2，…，第p个主成分代替原来的变量，其中p≤n，n为自然数；
C5）根据计算主成分载荷，其中lij为主成分载荷，为行为i列为j的单位矩阵；
C6）根据得到各指标的权重值：其中lij为主成分载荷，/>为各指标的权重值，p为选定的主成分个数；
步骤3）中所述的植被覆盖与管理因子的提取是基于植被覆盖度，通过对建立的模型进行修正后，按照此模型公式计算得到，具体是指用回归分析方法建立了C与植被覆盖度FVC以百分数表示之间的关系式，当FVC大于9.632%时，将其带入公式求得的C值不大于1，具体根据
得到植被覆盖与管理因子，其中C为植被覆盖与管理因子，FVC为植被覆盖度；
所述步骤3）中水土保持工程措施因子的获取，具体是指林地、草地、城镇、农村居民点以及未利用地般不采取任何保护措施，其P值为1；耕地确定其值的大小需要大量的观测小区实测资料，对采用措施前后的土壤流失量进行观测，分析观测数值得到不同水土保持措施下P值的大小。
CN202110722290.0A 2021-06-29 2021-06-29 基于多源数据和多时相数据的水土流失定量监测方法 Active CN113591572B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110722290.0A CN113591572B (zh) 2021-06-29 2021-06-29 基于多源数据和多时相数据的水土流失定量监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110722290.0A CN113591572B (zh) 2021-06-29 2021-06-29 基于多源数据和多时相数据的水土流失定量监测方法

Publications (2)

Publication Number Publication Date
CN113591572A CN113591572A (zh) 2021-11-02
CN113591572B true CN113591572B (zh) 2023-08-15

Family

ID=78244834

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110722290.0A Active CN113591572B (zh) 2021-06-29 2021-06-29 基于多源数据和多时相数据的水土流失定量监测方法

Country Status (1)

Country Link
CN (1) CN113591572B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114036777B (zh) * 2021-11-26 2022-10-18 成都理工大学 一种浅层土质滑坡的早期识别方法
CN116090858B (zh) * 2022-11-08 2023-08-01 北京师范大学 水资源及坡度双重限制下的生态修复潜力评价方法及系统
CN115993336B (zh) * 2023-03-23 2023-06-16 山东省水利科学研究院 一种输水渠渠道两侧植被损坏监测方法及预警方法
CN116523150B (zh) * 2023-07-05 2023-09-26 南昌工程学院 一种林下水土流失监测及防治模拟的方法
CN117151430B (zh) * 2023-10-30 2024-01-30 水利部交通运输部国家能源局南京水利科学研究院 一种小流域水土保持治理优先度遥感评估方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103940974A (zh) * 2014-02-19 2014-07-23 西北农林科技大学 基于gis的中尺度流域土壤侵蚀时空动态分析方法
CN105004725A (zh) * 2015-08-04 2015-10-28 珠江水利委员会珠江水利科学研究院 一种水土保持综合治理土壤侵蚀变化量实时定量监测方法
CN107145476A (zh) * 2017-05-23 2017-09-08 福建师范大学 一种基于改进tf‑idf关键词提取算法
CN109581412A (zh) * 2019-01-17 2019-04-05 合肥工业大学 一种快速进行水土流失动态变化遥感监测的方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106295576B (zh) * 2016-08-12 2017-12-12 中国水利水电科学研究院 一种基于自然地理特征的水源类型解析方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103940974A (zh) * 2014-02-19 2014-07-23 西北农林科技大学 基于gis的中尺度流域土壤侵蚀时空动态分析方法
CN105004725A (zh) * 2015-08-04 2015-10-28 珠江水利委员会珠江水利科学研究院 一种水土保持综合治理土壤侵蚀变化量实时定量监测方法
CN107145476A (zh) * 2017-05-23 2017-09-08 福建师范大学 一种基于改进tf‑idf关键词提取算法
CN109581412A (zh) * 2019-01-17 2019-04-05 合肥工业大学 一种快速进行水土流失动态变化遥感监测的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于Landsat影像的多时相植被覆盖度与地形因子关系研究――以平潭岛为例;刘尧文;沙晋明;;福建师范大学学报(自然科学版)(第04期);全文 *

Also Published As

Publication number Publication date
CN113591572A (zh) 2021-11-02

Similar Documents

Publication Publication Date Title
CN113591572B (zh) 基于多源数据和多时相数据的水土流失定量监测方法
Le et al. Adequacy of satellite-derived precipitation estimate for hydrological modeling in Vietnam basins
Boilley et al. Comparison between meteorological re-analyses from ERA-Interim and MERRA and measurements of daily solar irradiation at surface
Fang et al. Spatial downscaling of TRMM precipitation data based on the orographical effect and meteorological conditions in a mountainous area
Millward et al. Adapting the RUSLE to model soil erosion potential in a mountainous tropical watershed
Han et al. A soil moisture estimation framework based on the CART algorithm and its application in China
Wang et al. Development of a large-scale remote sensing ecological index in arid areas and its application in the Aral Sea Basin
Kim et al. Statistical downscaling for daily precipitation in Korea using combined PRISM, RCM, and quantile mapping: Part 1, methodology and evaluation in historical simulation
Joshi et al. Rainfall-Runoff Simulation in Cache River Basin, Illinois, Using HEC-HMS
Desalegn et al. Developing GIS-based soil erosion map using RUSLE of Andit Tid Watershed, central highlands of Ethiopia
CN114201922A (zh) 基于InSAR技术的动态滑坡敏感性预测方法及系统
Baik et al. Agricultural drought assessment based on multiple soil moisture products
Salmani-Dehaghi et al. Spatiotemporal assessment of the PERSIANN family of satellite precipitation data over Fars Province, Iran
Pachri et al. Development of water management modeling by using GIS in Chirchik river basin, Uzbekistan
Fitzgerald et al. Directed sampling using remote sensing with a response surface sampling design for site-specific agriculture
Wang et al. Evaluation of multi‐source precipitation data in a watershed with complex topography based on distributed hydrological modeling
Leng et al. Determination of all-sky surface soil moisture at fine spatial resolution synergistically using optical/thermal infrared and microwave measurements
CN117494419A (zh) 一种多模型耦合的流域水土流失遥感监测方法
Jiang et al. Characterizing precipitation uncertainties in a high-altitudinal permafrost watershed of the Tibetan plateau based on regional water balance and hydrological model simulations
Chantraket et al. An operational weather radar-based calibration of Z–R relationship over central region of Thailand
Rossi et al. Potential of MERIS fAPAR for drought detection
Li et al. Improvement of the multi-source weighted-ensemble precipitation dataset and application in the arid area of Tianshan Mountains, central Asia
Yang et al. Seasonal prediction of crop yields in Ethiopia using an analog approach
Sui et al. Disentangling error structures of precipitation datasets using decision trees
Gado et al. Evaluation of satellite-based rainfall estimates in the upper Blue Nile basin

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant