CN115358507A - 生产建设项目扰动图斑水土流失风险识别评估方法 - Google Patents

生产建设项目扰动图斑水土流失风险识别评估方法 Download PDF

Info

Publication number
CN115358507A
CN115358507A CN202210765136.6A CN202210765136A CN115358507A CN 115358507 A CN115358507 A CN 115358507A CN 202210765136 A CN202210765136 A CN 202210765136A CN 115358507 A CN115358507 A CN 115358507A
Authority
CN
China
Prior art keywords
data
pattern spot
risk
disturbance pattern
area
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.)
Granted
Application number
CN202210765136.6A
Other languages
English (en)
Other versions
CN115358507B (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.)
Water And Oil Maintenance Monitoring Center Ministry Of Water Resources
Pearl Water Soil And Water Conservation Monitoring Station Pearl Water Resources Commission
Pearl River Hydraulic Research Institute of PRWRC
Original Assignee
Water And Oil Maintenance Monitoring Center Ministry Of Water Resources
Pearl Water Soil And Water Conservation Monitoring Station Pearl Water Resources Commission
Pearl River Hydraulic Research Institute of PRWRC
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 Water And Oil Maintenance Monitoring Center Ministry Of Water Resources, Pearl Water Soil And Water Conservation Monitoring Station Pearl Water Resources Commission, Pearl River Hydraulic Research Institute of PRWRC filed Critical Water And Oil Maintenance Monitoring Center Ministry Of Water Resources
Priority to CN202210765136.6A priority Critical patent/CN115358507B/zh
Publication of CN115358507A publication Critical patent/CN115358507A/zh
Application granted granted Critical
Publication of CN115358507B publication Critical patent/CN115358507B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/10—Services
    • G06Q50/26—Government or public services
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/26—Segmentation of patterns in the image field; Cutting or merging of image elements to establish the pattern region, e.g. clustering-based techniques; Detection of occlusion
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/20—Image preprocessing
    • G06V10/28—Quantising the image, e.g. histogram thresholding for discrimination between background and foreground patterns
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V20/00—Scenes; Scene-specific elements
    • G06V20/10—Terrestrial scenes

Landscapes

  • Engineering & Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • Theoretical Computer Science (AREA)
  • Economics (AREA)
  • Multimedia (AREA)
  • Tourism & Hospitality (AREA)
  • Strategic Management (AREA)
  • Marketing (AREA)
  • Development Economics (AREA)
  • Educational Administration (AREA)
  • General Business, Economics & Management (AREA)
  • Entrepreneurship & Innovation (AREA)
  • General Health & Medical Sciences (AREA)
  • Primary Health Care (AREA)
  • Health & Medical Sciences (AREA)
  • Game Theory and Decision Science (AREA)
  • Operations Research (AREA)
  • Quality & Reliability (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了生产建设项目扰动图斑水土流失风险识别评估方法，收集研究区域生产建设项目扰动图斑矢量数据、坡度或数值高程模型DEM栅格数据以及高分辨遥感影像数据；根据以上数据计算水土流失风险识别9类评价指标VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR以及VMJ，从而获得水土流失风险识别值并判定风险等级。本发明能够及时发现问题、预测扰动图斑风险点的位置及其风险等级，为水行政部门履行监管职责提供必要技术手段，提高监管效能。

Description

生产建设项目扰动图斑水土流失风险识别评估方法
技术领域
本发明属于水土保持、智慧水利技术领域，涉及生产建设项目扰动图斑水土流失风险识别评估方法。
背景技术
我国连续3年组织开展水土保持遥感信息化监管工作，通过人机交互的方式提取扰动地块和合规性判读，并开展疑似违规图斑的外业复核、认定和查处，在生产建设项目水土保持遥感监管工作中，逐渐实现了遥感影像扰动地块解译技术从试点到深入应用。但人为水土流失类型多样，情况复杂，使得人为水土流失监管任务十分艰巨，国内对人为水土流失地块的自动化提取研究处于起步阶段，尤其是生产建设项目造成的水土流失监管风险预警模型尚未建立。对标部级遥感监管要求和智慧水土保持建设需求，目前主要有两方面技术问题亟待解决。一是要解决水土流失扰动地块人工提取与判别精度、效率不高的问题。目前遥感解译工作主要依托人机交互的方式开展，无法实现自动化智能解译，解译效率低、工作成本高。二是要解决人为水土流失风险难以发现和预警的问题。当前监管工作难以准确识别人为水土流失风险，需要有专业的模型支撑风险识别与预警，在生态环境监管日益严格的背景下，对于人为水土流失风险预警模型的研发和推广应用非常迫切。
发明内容
本发明的目的是提供一种生产建设项目扰动图斑水土流失风险识别评估方法，致力解决生产建设项目扰动图斑水土流失风险识别评估标准不统一、工作效率低下、投入成本较高等问题相关问题。
本发明所采用的技术方案是，生产建设项目扰动图斑水土流失风险识别评估方法，具体操作步骤如下：
步骤1：收集研究区域生产建设项目扰动图斑矢量数据；收集坡度或数值高程模型DEM栅格数据，对上述栅格数据进行重采样，采样后分辨率优于5m；采集研究区高分辨遥感栅格数据，遥感影像数据分辨率优于2m，包括RGB三波段；
步骤2：按扰动图斑矢量边界裁剪坡度或数值高程模型DEM栅格数据，裁剪结果记为数据A，根据裁剪后的数据A计算研究区域的扰动图斑内部坡度或高程均值VSmean、扰动图斑内部坡度或高程方差VSvar；
步骤3：沿原扰动图斑区域边界依次向外缓冲15米、30米和45米，分别得到缓冲区A、缓冲区B和缓冲区C；计算扰动图斑内部到外部坡度或高程连续降低点位比例VDecR、扰动图斑内外部坡度或高程急剧变化点位占比VSteR、扰动图斑内外部坡度或高程下降比例均值VDecM；
步骤4：按照扰动图斑区域边界向外缓冲500m，得到缓冲区D，使用缓冲区D坡度或数值高程模型DEM栅格数据，裁剪结果记为数据B,根据裁剪后的数据B计算扰动图斑所在区域坡度或高程均值VLSM、扰动图斑所在区域坡度或高程方差VLSVa；
步骤5：根据步骤1获得的研究区高分辨遥感影像数据，采用波段运算得到研究区单波段遥感影像灰度图，根据单波段遥感影像灰度图像元值计算扰动图斑裸露面积占比VBareR，在步骤1的扰动图斑矢量数据中提取扰动图斑正射投影面积VMJ；
步骤6：对步骤2-5得到的扰动图斑水土流失风险识别9类评价指标VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ对进行重要性排序，采用层次分析法或灰色隶属度判别法，定量计算各风险识别评价指标权重因子，分别为：wSmean、wSvar、wDecR、wSteR、wDecM、wLSM、wLSVa、wBaredD、wMJ；根据上述评价指标以及各风险识别评价指标权重因子确定水土流失风险识别值，风险识别值越大水土流失风险越大；具体风险等级判定标准如下：
风险识别值为
Figure BDA0003724523620000031
基本无风险；
风险识别值为
Figure BDA0003724523620000032
低风险；
风险识别值为
Figure BDA0003724523620000033
中风险；
风险识别值为
Figure BDA0003724523620000034
高风险；
风险识别值为
Figure BDA0003724523620000035
极高风险。
本发明的特点还在于，
步骤2具体如下：
按扰动图斑矢量边界裁剪坡度或数值高程模型DEM栅格数据，裁剪结果记为数据A，所述数据A包括M个像元，ki即为该像元对应的坡度或高程值，则有：
VSmean＝(k1+k2+k3+…+kN)/M (1)
VSvar＝[(k1-VSmean)2+(k2-VSmean)2+(k3-VSmean)2+…+(kN-VSmean)2+]/M (2)。
步骤3具体如下：
对原扰动图斑区域依次向外缓冲15米、30米和45米，分别得到缓冲区A、缓冲区B和缓冲区C，将缓冲区C最大轮廓边界正交划分为4个区域，即生成2×2正交单元格矢量数据，4个区域依次为区域1、区域2、区域3、区域4；
区域1、区域2、区域3、区域4分别与扰动图斑、缓冲区A、缓冲区B和缓冲区C进行空间交集运算由此共得到16个分区；
将每个分区对应的坡度或数值高程值，记为Gi-j，其中i表示区域1-4，分别用1、2、3、4表示；j表示扰动图斑、缓冲区A、缓冲区B、缓冲区C，分别用1、2、3、4表示；
1)令X1＝0，X1用于表示区域i的坡度或海拔高程连续降低的个数；即Gi-1、Gi-2、Gi-3、Gi-4连续降低的个数；
如果Gi-1>Gi-2>Gi-3>Gi-4，表示第i个区域的Gi-1、Gi-2、Gi-3、Gi-4是连续降低的，水土流失风险巨大，则X1累加1；同时计算提取“分区i-1”的中心坐标经纬度为Riskpoint(xi，yi)，该坐标点也即为该扰动图斑水土流失风险点；
则扰动图斑内部到外部坡度或高程连续降低点位比例VDecR＝X1/4 (3)
2)令X2＝0，X2用于表示区域i的坡度或海拔高程急剧降低的个数，
如果(Gi-1-Gi-4)/Gi-1>0.5，表示区域i坡度或海拔高程剧烈变化比例超过50％，水土流失风险巨大，则X2累加1；
则扰动图斑内外部坡度或高程急剧变化点位占比VSteR＝X2/4 (4)
扰动图斑内外部坡度或高程下降比例均值DecM的计算方法如下：
R1＝((G1-1-G1-2)/G1-1+(G1-2-G1-3)/G1-2+(G1-3-G1-4)/G1-3)/3 (6)
R2＝((G2-1-G2-2)/G2-1+(G2-2-G2-3)/G2-2+(G2-3-G2-4)/G2-3)/3 (7)
R3＝((G3-1-G3-2)/G3-1+(G3-2-G3-3)/G3-2+(G3-3-G3-4)/G3-3)/3 (8)
R4＝((G4-1-G4-2)/G4-1+(G4-2-G4-3)/G4-2+(G4-3-G4-4)/G4-3)/3 (9)
Ri表示第“区域i”坡度连续变化比率的均值，i＝1、2、3、4。
VDecM＝(R1+R2+R3+R4)/4 (10)。
扰动图斑所在大区域坡度/高程均值VLSM、扰动图斑所在大区域坡度/高程方差VLSVa的计算方法为：
按照扰动图斑边界向外缓冲500m，得到缓冲区D，使用缓冲区D裁剪坡度或数值高程模型DEM栅格数据，裁剪结果记为数据C，缓冲区D坡度或高程均值及方差越大，扰动图斑水土流失风险越高；
其中数据C包括U个像元，每个像元值hi即为该像元对应的坡度或高程值，则有：
VLSM＝(h1+h2+h3+…+hM)/U (11)
VLSVa＝[(h1-VLSM)2+(h2-VLSM)2+(h3-VLSM)2+…+(hU-VLSM)2]/U (12)。
步骤5具体如下：
使用扰动图斑裁剪研究区高分辨遥感影像栅格数据，裁剪结果记为数据D，其中数据D包括RGB三个波段，使用bandr、bandg和bandb分别表示数据D三个波段各像元值；
首先采用如下公式得到数据D的单波段灰度图记为bandgray，该数据记为数据E，
bandgray＝2×bandr-bandg-bandb (13)
其中，数据E各像元值范围为0到255。
设vthreshold为数据E二值化处理的临界阈值，T为数据E的像元个数，vi表示数据E第i个像元值；
当vi≥vthreshold时，表示第i个像元为非裸露，这种像元个数为T1；
当vi<vthreshold时，表示第i个像元为裸露，这种像元个数为T2；
显然T1+T2＝T；
则VBareR的计算公式为：
VBareR＝T2/T (14)。
对VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ进行归一化处理到[0,1]区间，分别得到归一化指标值为NORSmean、NORSvar、NORDecR、NORSteR、NORDecM、NORLSM、NORLSVa、NORBareR、NORMJ，那么水土流失风险识别值RiskV计算公式如下：
RiskV＝wSmean×NORSmean+wSvar×NORSvar+wDecR×NORDecR+wSteR×NORSteR+wDecM×NORDecM+wLSM×NORLSM+wLSVa×NORLSVa+wBareR×NORBareR+wMJ×NORMJ (15)
wSmean、wSvar、wDecR、wSteR、wDecM、wLSM、wLSVa、wBareR、wMJ分别为各评价指标的权重，权重值范围均为[0，1]，且有：wSmean+wSvar+wDecR+wSteR+wDecM+wLSM+wLSVa+wBareR+wMJ＝1；
显然，RiskV值越大扰动图斑水土流失风险越大，且RiskV值的变化范围为[0，1]。
每个分区对应的坡度或数值高程模型像元均值Gi-j的计算方法如下：
分别用16个分区裁剪步骤1的研究区坡度或数值高程模型DEM栅格数据，裁剪结果分别记为数据Bi-j，数据Bi-j包含Ni-j个像元，每个像元值px即为该像元对应的坡度或海拔高程值，x∈{1，2，3…Ni-j}；则数据Bi-j的坡度或海拔高程值Gi-j为：
Gi-j＝(p1+p2+p3+…+pNi-j)/Ni-j (3)
其中，Bi-j、Ni-j和Gi-j中的下标i＝1、2、3、4，分别表示区域1、区域2、区域3和区域4；下标j＝1、2、3、4，分别表示扰动图斑区域、缓冲区A、缓冲区B和缓冲区C。
对VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ进行归一化处理分别采用如下方法：
xi＝(Xi-Xmin)/(Xmax-Xmin)
其中，Xi为归一化处理之前原始数据序列第i个数据的值，Xmin和Xmax为归一化处理之前原始数据序列最小值和最大值，xi为归一化处理之后数据序列第i个数据的值。
步骤1的扰动图斑矢量数据为Shapefile格式，坐标系为CGCS2000。
本发明的有益效果是：本发明能够及时发现问题、预测扰动图斑风险点的位置及其风险等级，为生产建设项目水土保持信息化监管提供技术支撑，提升信息化监管效能，具体如下：
(1)本发明可为生产建设项目水土保持信息化监管提供技术支撑。
(2)能够准确识别评估生产建设项目扰动图斑水土流失风险等级、确定水土流失风险点地理坐标，对于降低生产建设项目水土流失风险、水土流失风险点针对性防控具有重要意义。
本发明方法可为生产建设项目水土保持信息化监管提供技术手段、提高信息化监管效能，对于水行政主管部门更好的履行水土保持信息化监管具有十分重要的意义。
附图说明
图1是本发明生产建设项目扰动图斑水土流失风险识别评估方法的流程图；
图2是本发明的扰动图斑缓冲形成的缓冲区A、B、C的结构示意图；
图3是本发明的分区结构示意图；
图4是本发明的缓冲区D结构示意图；
图5是本发明实施例中扰动图斑结构示意图；
图6是本发明实施例中扰动图斑所在区域坡度栅格数据结构示意图；
图7是本发明实施例中扰动图斑所在区域高分辨率遥感影像栅格数据结构示意图；
图8是本发明实施例中分别使用扰动图斑1、2 3裁剪的数据A对应的3个栅格数据；
图9是本发明实施例中扰动图斑1的3个缓冲区的分区图；
图10是本发明实施例中扰动图斑2的3个缓冲区的分区图；
图11是本发明实施例中扰动图斑3的3个缓冲区的分区图；
图12是本发明实施例中以扰动图斑1为边界向外缓冲得到的缓冲区D示意；
图13是本发明实施例中以扰动图斑1矢量裁剪的数据B示意图；
图14是图13的单波段灰度图；
图15是图14的像元值更新后二值化栅格数据图；
具体实施方式
下面结合附图对本发明的实施方式作出详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述实施例。
本发明的生产建设项目扰动图斑水土流失风险识别评估方法，其工作流程如图1所示，主要包括以下几个步骤：
一、扰动图斑水土流失风险识别评估指标确定
扰动图斑所在区域坡度或地形越复杂、扰动图斑裸露面积占比越大，该扰动图斑水土流失风险等级越高。基于地形地貌特征、扰动图斑裸露面积、扰动图斑面积等因素，确定了生产建设项目扰动图斑水土流失风险识别评估9类评价指标，具体分别为：VSmean(扰动图斑内部坡度或海拔高程均值)、VSvar(扰动图斑内部坡度或海拔高程方差)、VDecM(扰动图斑内部到外部坡度或海拔高程连续降低点位占比)、VDecR(内部到外部坡度或海拔高程急剧降低点位占比)、VDecM(扰动图斑内部到外部坡度或海拔高程下降比例均值)、VLSM(扰动图斑所在区域坡度或海拔高程均值)、VLSVa(扰动图斑所在区域坡度或海拔高程方差)、VBareR(扰动图斑裸露面积占扰动图斑面积比例)、VMJ(扰动图斑面积)。
上述9类评价指标均为正向指标，即指标值越大扰动图斑水土流失风险越大。二、扰动图斑水土流失风险识别评估指标权重确定
根据上述9类评价指标对扰动图斑水土流失风险影响的重要性进行由大到小排序，采用层次分析法或灰色隶属度判别法计算上述9类评价指标权重值大小，分别为：wSmean、wSvar、wDecR、wSteR、wDecM、wLSM、wLSVa、wBareR、wMJ。各指标权重值大小范围均为为[0，1]，且有：wSmean+wSvar+wDecR+wSteR+wDecM+wLSM+wLSVa+wBareR+wMJ＝1。
三、基础数据准备与处理
1.研究区域生产建设项目扰动图斑矢量数据(Shapefile格式，坐标系为CGCS2000)，矢量数据属性字段无特别要求；
2.研究区坡度或数值高程模型(DEM)栅格数据(GeoTiff格式，分辨率优于30m)，该栅格数据仅有1个波段，栅格数据各像元值为坡度或海拔高程。为保证计算精度等要求，需要对上述栅格数据进行重采样，采样后分辨率优于5m。
3.研究区高分辨遥感影像栅格数据(GeoTiff格式，分辨率优于2m，RGB三波段)。
四、水土流失风险识别评估因子值计算
(1)VSmean、VSvar指标值计算：
按扰动图斑矢量边界裁剪“坡度或数值高程模型(DEM)栅格数据”，裁剪结果记为“数据A”。
“数据A”包括M个像元，每个像元值ki即为该像元对应的坡度或高程值，则有：
VSmean＝(k1+k2+k3+…+kN)/M (1)
VSvar＝[(k1-VSmean)2+(k2-VSmean)2+(k3-VSmean)2+…+(kN-VSmean)2+]/M (2)
(2)VDecR、VSteR、VDecM指标值计算：
对原扰动图斑区域依次向外缓冲15米、30米和45米，如下图的缓冲区A、缓冲区B和缓冲区C，如图2所示。
将缓冲区C最大轮廓边界正交划分为4个区域，即生成2×2正交单元格矢量数据，4个区域依次为区域1、区域2、区域3、区域4；区域1、区域2、区域3、区域4分别与扰动图斑、缓冲区A、缓冲区B和缓冲区C进行空间交集运算，分别得到1-1、1-2、1-3、1-4、2-1、2-2、2-3、2-4、3-1、3-2、3-3、3-4、4-1、4-2、4-3、4-4共计16个分区，如图3所示。
分别用1-1、1-2、1-3、1-4、2-1、2-2、2-3、2-4、3-1、3-2、3-3、3-4、4-1、4-2、4-3、4-4共计16个分区裁剪研究区坡度或数值高程模型(DEM)栅格数据，裁剪结果分别记为“数据B1-1”、“数据B1-2”、“数据B1-3”、“数据B1-4”、“数据B2-1”、“数据B2-2”、“数据B2-3”、“数据B2-4”、“数据B3-1”、“数据B3-2”、“数据B3-3”、“数据B3-4”、“数据B4-1”、“数据B4-2”、“数据B4-3”、“数据B4-4”。
“数据B1-1”包含N1-1个像元，每个像元值pi即为该像元对应的坡度或海拔高程值，则“数据B1-1”的坡度或海拔高程值G1-1(单位：度或米)：
G1-1＝(p1+p2+p3+…+pN)/N1-1
同理计算“数据B1-2”、“数据B1-3”、“数据B1-4”、“数据B2-1”、“数据B2-2”、“数据B2-3”、“数据B2-4”、“数据B3-1”、“数据B3-2”、“数据B3-3”、“数据B3-4”、“数据B4-1”、“数据B4-2”、“数据B4-3”、“数据B4-4”数据的坡度或海拔高程值G1-2、G1-3、G1-4、G2-1、G2-2、G2-3、G2-4、G3-1、G3-2、G3-3、G3-4、G4-1、G4-2、G4-3、G4-4。
上述“数据Bi-j”、Ni-j和Gi-j中的下标i＝1、2、3、4，分别表示“区域1”、“区域2”、“区域3”和“区域4”；下标j＝1,、2、3、4，分别表示“扰动图斑区域”、“缓冲区A”、“缓冲区B”和“缓冲区C”。
那么，VDecR、VSteR、VDecM指标值计算流程如下：
1)令X1＝0，X1用于表示“区域i”的坡度或海拔高程连续降低的个数，即为“分区i-1”、“分区i-2”、“分区i-3”、“分区i-4”的Gi-1、Gi-2、Gi-3、Gi-4连续降低的个数。
如果Gi-1>Gi-2>Gi-3>Gi-4，则表示“区域i”的“分区i-1”、“分区i-2”、“分区i-3”、“分区i-4”的Gi-1、Gi-2、Gi-3、Gi-4是连续降低的，该分区在面临自然降雨后容易产生水土流失，且水土流失对扰动图斑外面有重要影响作用，该扰动图斑该分区水土流失风险巨大，让X1累加1。同时计算提取“分区i-1”的中心坐标经纬度为Riskpoint(xi，yi)，该坐标点也即为该扰动图斑水土流失风险点。
则VDecR指标值计算公式为：
VDecR＝X1/4 (4)
2)令X2＝0，X2用于表示“区域i”的坡度或海拔高程急剧降低的个数。
如果(Gi-1-Gi-4)/Gi-1>0.5，表示“区域i”坡度或海拔高程剧烈变化比例超过50％，水土流失风险巨大，则X2累加1；
则VSteR指标值计算公式为：
VSteR＝X2/4 (5)
3)VDecM表示扰动图斑与外部缓冲区A、B、C坡度或海拔高程变化比例的均值，具体计算公式如下：
R1＝((G1-1-G1-2)/G1-1+(G1-2-G1-3)/G1-2+(G1-3-G1-4)/G1-3)/3 (6)
R2＝((G2-1-G2-2)/G2-1+(G2-2-G2-3)/G2-2+(G2-3-G2-4)/G2-3)/3 (7)
R3＝((G3-1-G3-2)/G3-1+(G3-2-G3-3)/G3-2+(G3-3-G3-4)/G3-3)/3 (8)
R4＝((G4-1-G4-2)/G4-1+(G4-2-G4-3)/G4-2+(G4-3-G4-4)/G4-3)/3 (9)
Ri表示第“区域i”坡度连续变化比率的均值，i＝1、2、3、4。
则VDecM指标值计算公式为：
VDecM＝(R1+R2+R3+R4)/4 (10)。
(3)VLSM、VLSVa指标值计算：
按照扰动图斑边界向外缓冲500m，得到“缓冲区D”(如图5所示)。使用“缓冲区D”裁剪“坡度或数值高程模型(DEM)栅格数据”，裁剪结果记为“数据C”。显然，“缓冲区D”坡度或海拔高程均值及方差越大，表示扰动图斑所处区域大地形越复杂，该扰动图斑更容易发生水土流失风险，扰动图斑水土流失风险也越高。
“数据C”包括U个像元，每个像元值hi即为该像元对应的坡度或高程值，则VLSM、VLSVa指标值计算公式为：
VLSM＝(h1+h2+h3+…+hM)/U (11)
VLSVa＝[(h1-VLSM)2+(h2-VLSM)2+(h3-VLSM)2+…+(hU-VLSM)2]/U (12)
(4)VBareD指标值计算：
VBareD指标表示扰动图斑裸露面积占图斑面积的比例，显然该值越大扰动图斑水土流失风险越大。
使用扰动图斑裁剪“高分辨遥感影像栅格数据”，裁剪结果记为“数据D”，显然“数据D”，包括RGB三个波段，使用bandr、bandg和bandb分别表示“数据D”三个波段各像元值。
首先采用如下公式得到“数据D”的单波段灰度图(记为bandgray)，该数据记为“数据E”。
bandgray＝2×bandr-bandg-bandb (13)
其中“数据E”各像元值范围为0到255。
设vthreshold为“数据E”二值化处理的临界阈值，T为“数据E”的像元个数，vi表示“数据E”第i个像元值。
当vi≥vthreshold时，表示第i个像元为非裸露，这种像元个数为T1；
当vi<vthreshold时，表示第i个像元为裸露，这种像元个数为T2；
显然T1+T2＝T。
则VBareD指标的计算公式为：
VBareD＝T2/T (14)
(5)VMJ指标值计算：
扰动图斑正射投影面积即为VMJ。
五、水土流失风险识别计算
采用归一化公式对VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ进行归一化处理，分别得到归一化指标值为NORSmean、NORSvar、NORDecR、NORSteR、NORDecM、NORLSM、NORLSVa、NORBareR、NORMJ，那么水土流失风险识别值RiskV计算公式如下：
RiskV＝wSmean×NORSmean+wSvar×NORSvar+wDecR×NORDecR+wSteR×NORSteR+wDecM×NORDecM+wLSM×NORLSM+wLSVa×NORLSVa+wBareR×NORBareR+wMJ×NORMJ(15)
显然，RiskV值越大扰动图斑水土流失风险越大，且RiskV值的变化范围为[0，1]。
六、水土流失风险评估计算
RiskV参考值以及风险等级参照表1风险等级判定标准。
表1生产建设项目扰动图斑水土流失风险等级判定标准
Figure BDA0003724523620000131
Figure BDA0003724523620000141
具体实施例1
下面结合附图对本发明的具体实施方式作出详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述实施例。
一、扰动图斑水土流失风险识别评估指标确定
二、扰动图斑水土流失风险识别评估指标权重确定
根据上述9类评价指标对扰动图斑水土流失风险影像的重要性进行由大到小排序，重要性从大到小依次排序结果为：
VBareR>VSvar>VSmean>VDecR>VSteR>VDecM>VMJ>VLSVa>VLSM。
本具体实施例中采用层次分析法依次确定上述9个指标权重值大小分别为：wBareR＝0.2519、wSvar＝0.1385、wSmean＝0.1133、wDecR＝0.1053、wSteR＝0.1026、wDecM＝0.0887、wMJ＝0.0831、wLSVa＝0.0641、wLSM＝0.0525。
各指标权重值大小范围均为为[0，1]，且有：wSmean+wSvar+wDecR+wSteR+wDecM+wLSM+wLSVa+wBareR+wMJ＝1。
三、基础数据准备、预处理
1.某生产建设项目扰动图斑矢量数据(Shapefile格式，坐标系为CGCS2000)，图斑总个数为3个，具体如图5所示。
2.扰动图斑所在区域坡度栅格数据(记为“数据A”)，如图6所示；扰动图斑所在区域高分辨率遥感影像栅格数据(记为“数据B”)，如图7所示。其中坡度栅格数据分辨率为5m，满足计算要求。高分辨率遥感影像栅格数据分辨率为2m，满足计算要求。
四、水土流失风险识别评估因子值计算
(1)VSmean、VSvar指标值计算：
分别使用“扰动图斑1”、“扰动图斑2”、“扰动图斑3”裁剪“数据A”，分别得到3个栅格数据(GeoTiff格式)，依次记为“数据A1”、“数据A2”、“数据A3”，如图8所示。
经计算，“数据A1”共有像元2410个，其像元最大值为58.81(度)、最小值为1.69(度)，所有像元平均值VSmean_A1＝21.44(度)、方差为VSvar_A1＝232.08(度2)。
“数据A2”共有像元1055个，其像元最大值为60.26(度)、最小值为5.43(度)，所有像元平均值VSmean_A1＝30.65(度)、方差为VSvar_A1＝275.59(度2)。
“数据A3”共有像元1806个，其像元最大值为25.69(度)、最小值为0.75(度)，所有像元平均值VSmean_A1＝13.04(度)、方差为VSvar_A1＝20.61(度2)。
因此，“扰动图斑1”、“扰动图斑2”、“扰动图斑3”Smean、Svar指标值分别如下表所示：
表1
V<sub>Smean</sub>/(度) V<sub>Svar</sub>/(度<sup>2</sup>)
扰动图斑1 21.44 232.08
扰动图斑2 30.65 275.59
扰动图斑3 13.04 20.61
(2)VDecR、VSteR、VDecM指标值计算：
对“扰动图斑1”依次向外缓冲15米、30米和45米，分别得到“缓冲区A”、“缓冲区B”和“缓冲区C”。按照“缓冲区C”最大轮廓边界制作2×2正交单元格矢量数据，记为“正交单元格数据”。
“扰动图斑1”、“缓冲区A”、“缓冲区B”、“缓冲区C”和“正交单元格数据”均为矢量数据(Shapefile格式)，对上述5个数据进行两两空间交集运算，得到1-1、1-2、1-3、1-4、2-1、2-2、2-3、2-4、3-1、3-2、3-3、3-4、4-1、4-2、4-3、4-4共计16个分区，如图9所示。
使用上述16个分区裁剪“数据A”，依次得到“数据A1-1”、“数据A1-2”、“数据A1-3”、“数据A1-4”、A2-1”、“数据A2-2”、“数据A2-3”、“数据A2-4”、A3-1”、“数据A3-2”、“数据A3-3”、“数据A3-4”、A4-1”、“数据A4-2”、“数据A4-3”、“数据A4-4”，共计16个栅格数据(GeoTiff格式)。分别计算16个分区的平均坡度值，记为Gij，计算结果如表2所示。
表2“扰动图斑1”对应的16个分区的平均坡度值
Figure BDA0003724523620000161
Figure BDA0003724523620000171
因为有：
G11<G12<G13<G14
G21>G22<G23<G24
G31<G32<G33<G34
G41<G42<G43<G44，
所以，4个区域中有3个区域的坡度均值是连续降低的，故而有：
“扰动图斑1”的VDecR＝3/4＝0.75。
因为计算得到：
(G11-G14)/G11＝0.1914<0.5
(G21-G24)/G21＝-0.0566<0.5
(G31-G34)/G31＝0.6046>0.5
(G41-G44)/G41＝0.5939>0.5，
所以4个区域中有2个区域的坡度急剧降低的，故而有：
“扰动图斑1”的VSteR＝2/4＝0.50。
计算得到：
R1＝((G1-1-G1-2)/G1-1+(G1-2-G1-3)/G1-2+(G1-3-G1-4)/G1-3)/3＝0.0669
R2＝((G2-1-G2-2)/G2-1+(G2-2-G2-3)/G2-2+(G2-3-G2-4)/G2-3)/3＝-0.0197
R3＝((G3-1-G3-2)/G3-1+(G3-2-G3-3)/G3-2+(G3-3-G3-4)/G3-3)/3＝0.2220
R4＝((G4-1-G4-2)/G4-1+(G4-2-G4-3)/G4-2+(G4-3-G4-4)/G4-3)/3＝0.2241
所以有：
“扰动图斑1”的VDecM＝[R1+R2+R3+R4]/4＝0.1233。
同理，制作“扰动图斑2”和“扰动图斑3”的三个缓冲区和2×2正交单元格矢量数据(如图10和图11所示)。
同样的方法计算得到“扰动图斑2”和“扰动图斑3”的VDecR、VSteR、VDecM指标值，如表3所示。
表3
项目 V<sub>DecR</sub> V<sub>SteR</sub> V<sub>DecM</sub>
扰动图斑2 0.50 0.25 0.1057
扰动图斑3 0.50 0.5 -0.0984
(3)VLSM、VLSVa指标值计算：
以“扰动图斑1”为边界向外缓冲500米得到缓冲区矢量边界，记为“缓冲区D”(如图12所示)，用“缓冲区D”裁剪“数据A”得到栅格数据“数据Alarge_1”(GeoTiff格式)。
经计算，“数据Alarge_1”共有像元52117个，其像元最大值为60.26(度)、最小值为0.68(度)，所有像元平均值VSmean_A1＝19.36(度)、
方差为VSvar_A1＝159.26(度2)。
同样的计算方法得到“扰动图斑2”和“扰动图斑3”的VLSM、VLSVa指标值如下表所示。
表4
V<sub>LSM</sub>/(度) V<sub>LSVa</sub>/(度<sup>2</sup>)
扰动图斑2 30.39 269.7
扰动图斑3 16.24 99.59
(4)VBareD指标值计算：
使用“扰动图斑1”矢量文件裁剪“数据B”，裁剪结果记为“数据B1”，其中“数据B1”包括RGB三个波段，分别为bandr、bandg和bandb，如图13所示。
使用公式bandgray＝2×bandr-bandg-bandb计算得到“数据B1”的单波段灰度图，计算结果记为“数据B2”，如图14所示。“数据B2”像元值变化范围为[-30,158]。
经计算，“数据B2”的覆盖像元与裸露像元分割阈值为55.3，将“数据B2”像元值不高于55.31的像元赋值为0，像元高于55.31的像元赋值为1，像元值更新后的数据为二值化栅格数据，记为“数据B3”，结果如图15所示(黑色表示覆盖区域、灰色表示裸露区域)。其中像元值为0的像元有13585个，像元值为1的像元有3236个，因此，“扰动图斑1”的BareD指标值为：
VBareD＝13585/(13585+3236)＝0.8081
使用同样的方法计算得到“扰动图斑2”和“扰动图斑3”的BareD指标值分别为0.1619和0.1007。
(5)VMJ指标值计算：
借助ArcGIS等地理信息工具，经计算“扰动图斑1”、“扰动图斑2”、“扰动图斑3”的MJ指标值分别为60195.03m2、45140.14m2和26385.84m2。
综上，“扰动图斑1”、“扰动图斑2”、“扰动图斑3”的9个指标值计算结果如表5所示。
表5
Figure BDA0003724523620000191
Figure BDA0003724523620000201
五、水土流失风险识别计算
采用如下归一化公式对表5中的数据进行归一化处理到[0,1]区间。
xi＝(Xi-Xmin)/(Xmax-Xmin)
其中Xi为归一化处理之前原始数据序列第i个数据的值，Xmin和Xmax为归一化处理之前原始数据序列最小值和最大值，xi为归一化处理之后数据序列第i个数据的值。
表5中的数据归一化处理结果如表6所示。
表6
Figure BDA0003724523620000202
Figure BDA0003724523620000211
wSmean、wSvar、wDecR、wSteR、wDecM、wLSM、wLSVa、wBareR、wMJ分别为9个指标的权重值大小，NORSmean、NORSvar、NORDecR、NORSteR、NORDecM、NORLSM、NORLSVa、NORBareR、NORMJ为9个指标值，采用如下公式分别计算“扰动图斑1”、“扰动图斑2”、“扰动图斑3”的风险值。
RiskV＝wSmean×NORSmean+wSvar×NORSvar+wDecR×NORDecR+wSteR×NORSteR+wDecM×NORDecM+wLSM×NORLSM+wLSVa×NORLSVa+wBareR×NORBareR+wMJ×NORMJ
“扰动图斑1”、“扰动图斑2”、“扰动图斑3”的风险值依次为0.8022、0.5180和0.1026。
六、水土流失风险评估计算
根据生产建设项目扰动图斑水土流失风险等级判定标准，“扰动图斑1”、“扰动图斑2”、“扰动图斑3”风险等级分别为“极高风险”、“中风险”和“基本无风险”。

Claims (9)

1.生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，具体操作步骤如下：
步骤1：收集研究区域生产建设项目扰动图斑矢量数据；收集坡度或数值高程模型DEM栅格数据，对上述栅格数据进行重采样，采样后分辨率优于5m；采集研究区高分辨遥感栅格数据，遥感影像数据分辨率优于2m，包括RGB三波段；
步骤2：按扰动图斑矢量边界裁剪坡度或数值高程模型DEM栅格数据，裁剪结果记为数据A，根据裁剪后的数据A计算研究区域的扰动图斑内部坡度或高程均值VSmean、扰动图斑内部坡度或高程方差VSvar；
步骤3：沿原扰动图斑区域边界依次向外缓冲15米、30米和45米，分别得到缓冲区A、缓冲区B和缓冲区C；计算扰动图斑内部到外部坡度或高程连续降低点位比例VDecR、扰动图斑内外部坡度或高程急剧变化点位占比VSteR、扰动图斑内外部坡度或高程下降比例均值VDecM；
步骤4：按照扰动图斑区域边界向外缓冲500m，得到缓冲区D，使用缓冲区D坡度或数值高程模型DEM栅格数据，裁剪结果记为数据B,根据裁剪后的数据B计算扰动图斑所在区域坡度或高程均值VLSM、扰动图斑所在区域坡度或高程方差VLSVa；
步骤5：根据步骤1获得的研究区高分辨遥感影像数据，采用波段运算得到研究区单波段遥感影像灰度图，根据单波段遥感影像灰度图像元值计算扰动图斑裸露面积占比VBareR，在步骤1的扰动图斑矢量数据中提取扰动图斑正射投影面积VMJ；
步骤6：对步骤2-5得到的扰动图斑水土流失风险识别9类评价指标VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ对进行重要性排序，采用层次分析法或灰色隶属度判别法，定量计算各风险识别评价指标权重因子，分别为：wSmean、wSvar、wDecR、wSteR、wDecM、wLSM、wLSVa、wBaredD、wMJ；根据上述评价指标以及各风险识别评价指标权重因子确定水土流失风险识别值，风险识别值越大水土流失风险越大；具体风险等级判定标准如下：
风险识别值为
Figure FDA0003724523610000021
基本无风险；
风险识别值为
Figure FDA0003724523610000022
低风险；
风险识别值为0.4～0.6，中风险；
风险识别值为0.6～0.8，高风险；
风险识别值为0.8～1，极高风险。
2.根据权利要求1所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，步骤2具体如下：
按扰动图斑矢量边界裁剪坡度或数值高程模型DEM栅格数据，裁剪结果记为数据A，所述数据A包括M个像元，ki即为该像元对应的坡度或高程值，则有：
VSmean＝(k1+k2+k3+…+kN)/M (1)
VSvar＝[(k1-VSmean)2+(k2-VSmean)2+(k3-VSmean)2+…+(kN-VSmean)2+]/M (2)。
3.根据权利要求1所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，步骤3具体如下：
对原扰动图斑区域依次向外缓冲15米、30米和45米，分别得到缓冲区A、缓冲区B和缓冲区C，将缓冲区C最大轮廓边界正交划分为4个区域，即生成2×2正交单元格矢量数据，4个区域依次为区域1、区域2、区域3、区域4；
区域1、区域2、区域3、区域4分别与扰动图斑、缓冲区A、缓冲区B和缓冲区C进行空间交集运算由此共得到16个分区；
将每个分区对应的坡度或数值高程值，记为Gi-j，其中i表示区域1-4，分别用1、2、3、4表示；j表示扰动图斑、缓冲区A、缓冲区B、缓冲区C，分别用1、2、3、4表示；
1)令X1＝0，X1用于表示区域i的坡度或海拔高程连续降低的个数；即Gi-1、Gi-2、Gi-3、Gi-4连续降低的个数；
如果Gi-1＞Gi-2＞Gi-3＞Gi-4，表示第i个区域的Gi-1、Gi-2、Gi-3、Gi-4是连续降低的，水土流失风险巨大，则X1累加1；同时计算提取“分区i-1”的中心坐标经纬度为Riskpoint(xi，yi)，该坐标点也即为该扰动图斑水土流失风险点；
则扰动图斑内部到外部坡度或高程连续降低点位比例VDecR＝X1/4 (3)；
2)令X2＝0，X2用于表示区域i的坡度或海拔高程急剧降低的个数，
如果(Gi-1-Gi-4)/Gi-1＞0.5，表示区域i坡度或海拔高程剧烈变化比例超过50％，水土流失风险巨大，则X2累加1；
则扰动图斑内外部坡度或高程急剧变化点位占比VSteR＝X2/4 (4)
扰动图斑内外部坡度或高程下降比例均值DecM的计算方法如下：
R1＝((G1-1-G1-2)/G1-1+(G1-2-G1-3)/G1-2+(G1-3-G1-4)/G1-3)/3 (6)
R2＝((G2-1-G2-2)/G2-1+(G2-2-G2-3)/G2-2+(G2-3-G2-4)/G2-3)/3 (7)
R3＝((G3-1-G3-2)/G3-1+(G3-2-G3-3)/G3-2+(G3-3-G3-4)/G3-3)/3 (8)
R4＝((G4-1-G4-2)/G4-1+(G4-2-G4-3)/G4-2+(G4-3-G4-4)/G4-3)/3 (9)
Ri表示第“区域i”坡度连续变化比率的均值，i＝1、2、3、4。
VDecM＝(R1+R2+R3+R4)/4 (10)。
4.根据权利要求1所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，扰动图斑所在大区域坡度/高程均值VLSM、扰动图斑所在大区域坡度/高程方差VLSVa的计算方法为：
按照扰动图斑边界向外缓冲500m，得到缓冲区D，使用缓冲区D裁剪坡度或数值高程模型DEM栅格数据，裁剪结果记为数据C，缓冲区D坡度或高程均值及方差越大，扰动图斑水土流失风险越高；
其中数据C包括U个像元，每个像元值hi即为该像元对应的坡度或高程值，则有：
VLSM＝(h1+h2+h3+…+hM)/U (11)
VLSVa＝[(h1-VLSM)2+(h2-VLSM)2+(h3-VLSM)2+…+(hU-VLSM)2]/U (12)。
5.根据权利要求1所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，步骤5具体如下：
使用扰动图斑裁剪研究区高分辨遥感影像栅格数据，裁剪结果记为数据D，其中数据D包括RGB三个波段，使用bandr、bandg和bandb分别表示数据D三个波段各像元值；
首先采用如下公式得到数据D的单波段灰度图记为bandgray，该数据记为数据E，
bandgray＝2×bandr-bandg-bandb (13)
其中，数据E各像元值范围为0到255。
设vthreshold为数据E二值化处理的临界阈值，T为数据E的像元个数，vi表示数据E第i个像元值；
当vi≥vthreshold时，表示第i个像元为非裸露，这种像元个数为T1；
当vi＜vthreshold时，表示第i个像元为裸露，这种像元个数为T2；
显然T1+T2＝T；
则VBareR的计算公式为：
VBareR＝T2/T (14)。
6.根据权利要求1所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，对VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ进行归一化处理到[0,1]区间，分别得到归一化指标值为NORSmean、NORSvar、NORDecR、NORSteR、NORDecM、NORLSM、NORLSVa、NORBareR、NORMJ，那么水土流失风险识别值RiskV计算公式如下：
RiskV＝wSmean×NORSmean+wSvar×NORSvar+wDecR×NORDecR+wSteR×NORSteR+wDecM×NORDecM+wLSM×NORLSM+wLSVa×NORLSVa+wBareR×NORBareR+wMJ×NORMJ(15)
wSmean、wSvar、wDecR、wSteR、wDecM、wLSM、wLSVa、wBareR、wMJ分别为各评价指标的权重，权重值范围均为[0，1]，且有：wSmean+wSvar+wDecR+wSteR+wDecM+wLSM+wLSVa+wBareR+wMJ＝1；
显然，RiskV值越大扰动图斑水土流失风险越大，且RiskV值的变化范围为[0，1]。
7.根据权利要求3所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，每个分区对应的坡度或数值高程模型像元均值Gi-j的计算方法如下：
分别用16个分区裁剪步骤1的研究区坡度或数值高程模型DEM栅格数据，裁剪结果分别记为数据Bi-j，数据Bi-j包含Ni-j个像元，每个像元值px即为该像元对应的坡度或海拔高程值，x∈{1，2，3…Ni-j}；则数据Bi-j的坡度或海拔高程值Gi-j为：
Gi-j＝(p1+p2+p3+…+pNi-j)/Ni-j (3)
其中，Bi-j、Ni-j和Gi-j中的下标i＝1、2、3、4，分别表示区域1、区域2、区域3和区域4；下标j＝1、2、3、4，分别表示扰动图斑区域、缓冲区A、缓冲区B和缓冲区C。
8.根据权利要求6所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，对VSmean、VSvar、VDecR、VSteR、VDecM、VLSM、VLSVa、VBareR、VMJ进行归一化处理分别采用如下方法：
xi＝(Xi-Xmin)/(Xmax-Xmin)
其中，Xi为归一化处理之前原始数据序列第i个数据的值，Xmin和Xmax为归一化处理之前原始数据序列最小值和最大值，xi为归一化处理之后数据序列第i个数据的值。
9.根据权利要求1所述的生产建设项目扰动图斑水土流失风险识别评估方法，其特征在于，步骤1所述扰动图斑矢量数据为Shapefile格式，坐标系为CGCS2000。
CN202210765136.6A 2022-06-30 2022-06-30 生产建设项目扰动图斑水土流失风险识别评估方法 Active CN115358507B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210765136.6A CN115358507B (zh) 2022-06-30 2022-06-30 生产建设项目扰动图斑水土流失风险识别评估方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210765136.6A CN115358507B (zh) 2022-06-30 2022-06-30 生产建设项目扰动图斑水土流失风险识别评估方法

Publications (2)

Publication Number Publication Date
CN115358507A true CN115358507A (zh) 2022-11-18
CN115358507B CN115358507B (zh) 2023-05-05

Family

ID=84030865

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210765136.6A Active CN115358507B (zh) 2022-06-30 2022-06-30 生产建设项目扰动图斑水土流失风险识别评估方法

Country Status (1)

Country Link
CN (1) CN115358507B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117113038A (zh) * 2023-10-25 2023-11-24 珠江水利委员会珠江水利科学研究院 城市水土流失黄泥水事件溯源方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130090975A1 (en) * 2011-10-11 2013-04-11 Pioneer Hi-Bred International, Inc. Method for Managing a Cellulosic Biomass Harvest
CN110415346A (zh) * 2019-07-10 2019-11-05 华中师范大学 利用面向对象的三维元胞自动机进行水土流失模拟的方法
CN111707620A (zh) * 2020-06-11 2020-09-25 中国电建集团华东勘测设计研究院有限公司 土地利用分类规则集、水土流失监测方法及系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130090975A1 (en) * 2011-10-11 2013-04-11 Pioneer Hi-Bred International, Inc. Method for Managing a Cellulosic Biomass Harvest
CN110415346A (zh) * 2019-07-10 2019-11-05 华中师范大学 利用面向对象的三维元胞自动机进行水土流失模拟的方法
CN111707620A (zh) * 2020-06-11 2020-09-25 中国电建集团华东勘测设计研究院有限公司 土地利用分类规则集、水土流失监测方法及系统

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117113038A (zh) * 2023-10-25 2023-11-24 珠江水利委员会珠江水利科学研究院 城市水土流失黄泥水事件溯源方法及系统
CN117113038B (zh) * 2023-10-25 2024-02-09 珠江水利委员会珠江水利科学研究院 城市水土流失黄泥水事件溯源方法及系统

Also Published As

Publication number Publication date
CN115358507B (zh) 2023-05-05

Similar Documents

Publication Publication Date Title
WO2021258758A1 (zh) 一种基于多因素的海岸线变化识别方法
CN102521624B (zh) 一种土地利用类型分类的方法和系统
CN107067781B (zh) 一种用于先进驾驶辅助系统应用的gis道路黑点地图生成方法
CN110134907B (zh) 一种降雨缺失数据填补方法、系统及电子设备
CN110119556B (zh) 一种区域水源涵养功能的时空演变分析方法
CN110956207B (zh) 一种光学遥感影像全要素变化检测方法
CN114611834B (zh) 一种基于多维特征分析的电力发电站选址评估规划方法
CN112966941A (zh) 一种基于交通事故大数据的事故黑点识别方法及系统
CN115358507B (zh) 生产建设项目扰动图斑水土流失风险识别评估方法
CN116486289A (zh) 一种多源数据和知识驱动下的燃气管道高后果区识别方法
CN113158899B (zh) 基于遥感夜光暗目标增强技术的村镇发展状态测度方法
CN112836590B (zh) 洪涝灾害监测方法、装置、电子设备及存储介质
CN112907113B (zh) 一种考虑空间相关性的植被变化成因识别方法
CN112232303A (zh) 一种基于高分遥感影像的草原道路信息提取方法
CN110765900A (zh) 一种基于dssd的自动检测违章建筑方法及系统
CN113450461B (zh) 一种排泥库土工布点云提取方法
CN115457386A (zh) 一种村庄用地信息化生成方法
CN111476434B (zh) 一种基于gis的土壤重金属分形维数空间变异分析方法
CN113496182B (zh) 基于遥感影像的道路提取方法及装置、存储介质及设备
CN114612800A (zh) 一种核算城市群建筑物质存量及时空变化的方法、系统
CN113780880A (zh) 一种基于空间代表性的pm2.5观测站点布局评价方法
CN117849907B (zh) 基于多源数据的气象灾害靶向预警方法及系统
FAN et al. Intelligent antenna attitude parameters measurement based on deep learning ssd model
CN110751398B (zh) 一种区域生态质量评价方法及装置
CN113194152B (zh) 一种基于手机信令数据的低活动强度建设用地识别方法

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