CN110415346B - 利用面向对象的三维元胞自动机进行水土流失模拟的方法 - Google Patents

利用面向对象的三维元胞自动机进行水土流失模拟的方法 Download PDF

Info

Publication number
CN110415346B
CN110415346B CN201910619240.2A CN201910619240A CN110415346B CN 110415346 B CN110415346 B CN 110415346B CN 201910619240 A CN201910619240 A CN 201910619240A CN 110415346 B CN110415346 B CN 110415346B
Authority
CN
China
Prior art keywords
area
grade
neighbor
change
water
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
CN201910619240.2A
Other languages
English (en)
Other versions
CN110415346A (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.)
Central China Normal University
Original Assignee
Central China 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 Central China Normal University filed Critical Central China Normal University
Priority to CN201910619240.2A priority Critical patent/CN110415346B/zh
Publication of CN110415346A publication Critical patent/CN110415346A/zh
Application granted granted Critical
Publication of CN110415346B publication Critical patent/CN110415346B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20—Finite element generation, e.g. wire-frame surface description, tesselation

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Remote Sensing (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Processing Or Creating Images (AREA)

Abstract

本发明公开了一种利用面向对象的三维元胞自动机进行水土流失模拟的方法，计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图；计算各等级图斑面积和；预测研究区域内发生水土流失的总面积并得到各等级发生水土流失的面积；遍历土壤侵蚀强度等级矢量图所有图斑，获取图斑空间联合概率；根据坡度、植被覆盖度、土地利用类型及水保措施计算三维决策因子；根据“地理学第一定律”和空间关系，定义面向对象的三维元胞自动机的变化规则；基于发生等级变化后的土壤侵蚀强度等级图，重新计算空间联合概率，选取发生水土流失的区域。本发明以“地理学第一定律”为基础，提出了一种大尺度下面向对象的三维元胞自动机，实现水土流失的三维动态模拟。

Description

利用面向对象的三维元胞自动机进行水土流失模拟的方法
技术领域
本发明属于地理信息科学与地理过程模拟技术领域，涉及一种利用面向对象的三维元胞自动机进行水土流失模拟的方法，具体涉及一种基于三因子模型(地形坡度、地表覆盖度和土地利用类型)和不确定性的面向对象三维元胞自动机实现水土流失动态模拟的方法。
背景技术
现有水土流失模拟***需考虑地形、地表覆盖、降雨等多种影响因素。在建立水土流失***的过程中，由于参数多样，运算复杂，不利于水土流失的模拟。同时，当前水土流失模拟***中元胞自动机多以栅格为单元且为二维空间分析, 忽略了高程因素的影响，缺少三维可视化显示，不利于进行水土流失变化趋势的模拟；仅基于栅格数据建立模型，不能充分考虑地理实体的空间特征，对全局空间考虑不足。
发明内容
为应对三维空间对象的动态模拟问题，本发明以地理学第一定律为基础，提供了一种基于三因子模型(地形坡度、地表覆盖度和土地利用类型)和不确定性的面向对象三维元胞自动机实现水土流失动态模拟的方法。
本发明所采用的技术方案是：1.一种利用面向对象的三维元胞自动机进行水土流失模拟的***，其特征在于，包括以下步骤：
步骤1：计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图；
步骤2：遍历土壤侵蚀强度等级矢量图中所有图斑，得到ID号为IDcenter(i)、等级为j的图斑i的面积Aj center(i)；记录等级为j的所有图斑的面积和为Aj rank；
Figure BDA0002125003080000011
其中，j＝1,2,3,…,N，N为等级总数；
步骤3：预测研究区域内发生水土流失的总面积Asum，其中研究区域总面积为Aarea；
步骤4：遍历土壤侵蚀强度等级矢量图所有图斑，获取图斑空间联合概率；
步骤5：将相同研究区域内的植被覆盖度图及坡度图归一化后转换为栅格图；记x,y为栅格图行列号，z(x,y)为根据DEM数据得到的(x,y)像元处的高程值；结合水保措施情况，计算(x,y,z(x,y))像元的三维决策因子Pchange，Pchange越高更易发生水土流失；
Figure BDA0002125003080000021
Figure BDA0002125003080000022
Figure BDA0002125003080000023
其中，NDVI为植被指数，NDVImax、NDVImin为最大及最小NDVI值，V为植被覆盖度，L对应土地利用类型，Sx,y,z(x,y)为坡度，dz/dx和dz/dy分别表示像元及其邻域确定的x方向变化率和y方向变化率；W为水保措施情况，W＝0表示无水保措施，W＝1表示设置水保措施；
步骤6：根据“地理学第一定律”和空间关系，定义面向对象三维元胞的变化规则；
步骤7：对研究区域进行水文分析；
步骤8：基于发生等级变化后的土壤侵蚀强度等级图，重新计算空间概率，选取发生水土流失的区域。
作为优选，步骤1中，利用三因子方法计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图，等级数目为5；其中，所述三因子为地形坡度、地表覆盖度和土地利用类型。
作为优选，步骤3中，根据先验概率P(A)预测研究区域内发生水土流失的总面积Asum，其中研究区域总面积为Aarea；
由贝叶斯原理对发生水土流失的面积进行分配得出各等级j发生水土流失区域的面积Aj happen；
Figure BDA0002125003080000024
Asum＝Aarea×P(A)
Figure BDA0002125003080000031
Figure BDA0002125003080000032
其中，Aj rank为等级为j的所有图斑的面积和，j＝1,2,3,…,5；
Figure BDA0002125003080000033
为先验发生水土流失的面积比例，即先验概率；
P(Bj)表示发生第j等级水土流失的比例；
Figure BDA0002125003080000034
为先验发生水土流失时第j等级发生水土流失的面积比例。
作为优选，步骤4的具体实现包括以下子步骤：
步骤4.1：计算图斑i对应的邻接j等级图斑占其所有邻接图斑的数量比例 PLocal(Aij)与面积比例PLocal(Bij)，根据“地理学第一定律”，得到图斑i对应等级j 的空间局部概率Pj Local(i)、局部影响权重Wij，得出图斑i的局部最大影响等级 Rneighbor(i)及局部最大影响权重Wneighbor(i)；同时记录图斑i邻接图斑中等级为 Rneighbor(i)且面积最大的图斑对应的ID号为IDneighbor(i)，表示图斑i的最大影响邻接图斑；
Figure BDA0002125003080000035
Figure BDA0002125003080000036
Figure BDA0002125003080000037
Figure BDA0002125003080000038
其中，Pj Local(i)为图斑i对应等级j的空间局部概率，j＝1,2,3,…,5；
Wneighbor(i)＝max(Wij)
其中，图斑i局部最大影响权重对应的等级记为Rneighbor(i)，j＝1,2,3,…,5；
步骤4.2：将等级为j的图斑i的邻接图斑的数量和面积与全局中所有j等级图斑的邻接图斑数量和面积的最大值进行比较，得到数量比例PGlobal(Ai)与面积比例PGlobal(Bi)，计算图斑i的空间全局概率Pj Global(i)；
Figure BDA0002125003080000041
Figure BDA0002125003080000042
Figure BDA0002125003080000043
其中，j＝1,2,3,…,5。
作为优选，步骤6中所述面向对象三维元胞的变化规则，包括：
1、设定面向对象的选取规则并选择发生土壤侵蚀强度等级变化的图斑；
(1)如果图斑i的Wneighbor(i)＝1，该图斑的等级将变为Rneighbor(i)，记录其 IDcenter(i)于ChangeAll列表中，表示图斑内部全部发生等级变化；
(2)如果图斑i的Aj center(i)<a且Wneighbor(i)>b,该图斑内部等级全部变为 Rneighbor(i)，记录其IDcenter(i)于ChangeAll列表中；
(3)根据Pj Global(i)选取图斑内部部分区域发生等级变化的图斑，采用面积置信水平，数量置信水平及设置阈值三种方法进行计算；
将等级为j的所有图斑依据Pj Global(i)从大到小进行排序选取，面积选取为选取图斑面积和达到面积Aj choose时停止选取，数量选取为选取j等级图斑的设定比例数量的图斑，阈值选取为当Pj Global(i)大于设定阈值时进行选取；记录选取的图斑的IDcenter(i)于ChangePortion列表中表示图斑内部部分区域发生等级变化，公式如下：
Figure BDA0002125003080000044
其中，规则中a，b，c均代表置信水平，j＝1,2,3,…,5；
2、设定三维元胞自动机变化规则并改变土壤侵蚀强度等级；
将土壤侵蚀强度等级矢量图分别以图斑ID及图斑等级Rank为属性值分别转换为ID栅格及等级Rank栅格，将土地利用类型图转换为栅格；
遍历栅格，记录当前中心元胞为(x,y,z(x,y))，其中x,y为行列号，z(x,y)为根据DEM数据得到的(x,y)处高程值；t时刻中心元胞的ID为IDx,y,z(x,y)(t)，与其所位于的图斑i的ID号IDcenter(i)对应相等，土地利用类型为Lx,y,z(x,y)(t)，元胞的等级状态为Rx,y,z(x,y)(t)；以3×3为邻域定义x,y方向上的邻居元胞；其中，设定所有图斑 i的初始改变状态status(i)为0；
三维元胞自动机变化规则如下：
(1)当元胞(x,y,z(x,y))的IDx,y,z(x,y)(t)位于ChangeAll列表中，且三维决策因子Pchange不为0，即无水保措施，将该元胞的t+1时刻的等级Rx,y,z(x,y)(t+1)改变为 Rneighbor(i)；
(2)当元胞(x,y,z(x,y))与邻居元胞ID不相等时，即IDx,y,z(x,y)(t)不等于(x,y-1)位置的IDx,y-1,z(x,y-1)(t)或(x,y+1)位置的IDx,y+1,z(x,y+1)(t)时，元胞位于图斑的边界线；若IDx,y,z(x,y)(t)位于ChangePortion列表中且其邻域中存在邻居元胞ID号为IDneighbor(i) 时，记录像元位置为(m,n)，并以此为起始点像元开始进行局部变化；
局部等级变化规则有如下4种方法：
①在邻域内每次随机选择一个邻居元胞，依次生成连续成片的像元并将元胞 t+1时刻等级更改为Rneighbor(i)；
②选取邻居元胞中Pchange最大值更改等级为Rneighbor(i)；
③每次随机生成一个方向，如果其Pchange大于置信水平d，选中该元胞发生等级变化；若无邻居元胞符合条件，随机抽中一方向继续进行变化；
④计算邻居元胞Pchange占邻域范围所有邻居元胞Pchange和的百分比，即概率，以此为随机抽取的选中概率进行选取并改变等级；
选取过程中保证像元ID与起始点ID相等，即在该图斑内部发生变化，直到改变像元的面积和达到随机面积Achange(i)，设置status(i)为1；返回并由位置 (x,y+1)处像元开始继续遍历，公式如下：
Figure BDA0002125003080000051
其中，e为(0,1)之间的随机数,j＝1,2,3,…,5。
作为优选，步骤7中，根据DEM数据，确定元胞(x,y,z(x,y))位置处的流向Dirx,y,z(x,y)及流量Accx,y,z(x,y)；
Figure BDA0002125003080000061
Accx,y,z(x,y)＝F(Dirx,y,z(x,y),weight)
其中，n＝1,2,4,8,16,32,64,128，代表8个邻居元胞；dn表示元胞间距离，当 n＝1,4,16,64时，dn＝1；当n＝2,8,32,128时，dn＝1.5；流量函数F中Dirx,y,z(x,y)和 weight代表流向栅格和权重矩阵。
作为优选，步骤8中所述选取发生水土流失的区域，选取规则为对于j等级的所有图斑，以局部最大影响等级与局部最大影响权重为双重条件进行降序排序，图斑i的Rneighbor(i)与Wneighbor(i)越高，越易发生水土流失；按所排顺序选取图斑，直到选取的图斑面积和达到步骤3中计算的Aj happen为止；对步骤7得到的水文分析结果中提取出发生水土流失的区域，完成研究区域内水土流失模拟过程；其中，Aj happen是由贝叶斯原理对发生水土流失的面积进行分配得出各等级j发生水土流失区域的面积。
相对于现有技术，本发明的有益效果是：将三维决策因子加入到三维元胞自动机的变化规则中，考虑了水土流失中高程因素的影响；以空间局部概率与空间全局概率作为三维元胞自动机中等级变化的决策依据，综合考虑了面向对象的区域空间信息、全局空间信息和元胞局部空间信息，符合“地理学第一定律”；三维元胞自动机局部等级变化规则考虑了地理过程的不确定性影响；依据贝叶斯原理预测未来发生水土流失区域的面积，可实现水土流失的三维动态模拟。
附图说明
图1是本发明实施例的流程图；
图2是本发明实施例中面向对象选取规则示意图；
图3是本发明实施例中三维元胞自动机示意图；
图4是本发明实施例中基于DEM三维元胞自动机示意图；
图5是本发明实施例中元胞自动机x,y方向上邻域示意图；
图6是本发明实施例中流向方向代码示意图；
图7是本发明实施例中三维元胞自动机变化规则示意图。
具体实施方式
为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。
请见图1，本发明提供的一种利用面向对象的三维元胞自动机进行水土流失模拟的***，包括以下步骤：
步骤1：计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图；
本实施例中，利用三因子(地形坡度、地表覆盖度和土地利用类型)计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图，等级数目为5；
步骤2：遍历土壤侵蚀强度等级矢量图中所有图斑，得到ID号为IDcenter(i)、等级为j的图斑i的面积Aj center(i)；记录等级为j的所有图斑的面积和为Aj rank；
其中，j＝1,2,3,…,N，N为等级总数，本实施例为5；
Figure BDA0002125003080000071
步骤3：预测研究区域内发生水土流失的总面积Asum，其中研究区域总面积为Aarea；
本实施例中，根据先验概率P(A)预测研究区域内发生水土流失的总面积Asum，其中研究区域总面积为Aarea；
由贝叶斯原理对发生水土流失的面积进行分配得出各等级j发生水土流失区域的面积Aj happen；
Figure BDA0002125003080000072
Asum＝Aarea×P(A)
Figure BDA0002125003080000073
Figure BDA0002125003080000074
其中，Aj rank为等级为j的所有图斑的面积和，j＝1,2,3,…,5；
Figure BDA0002125003080000075
为先验发生水土流失的面积比例，即先验概率；
P(Bj)表示发生第j等级水土流失的比例(概率)；
Figure BDA0002125003080000081
为先验发生水土流失时第j等级发生水土流失的面积比例(概率)。
步骤4：遍历土壤侵蚀强度等级矢量图所有图斑，获取图斑空间联合概率；
步骤4.1：计算图斑i对应的邻接j等级图斑占其所有邻接图斑的数量比例 PLocal(Aij)与面积比例PLocal(Bij)，根据“地理学第一定律”，得到图斑i对应等级j 的空间局部概率Pj Local(i)、局部影响权重Wij，得出图斑i的局部最大影响等级 Rneighbor(i)及局部最大影响权重Wneighbor(i)；同时记录图斑i邻接图斑中等级为 Rneighbor(i)且面积最大的图斑对应的ID号为IDneighbor(i)，表示图斑i的最大影响邻接图斑；
Figure BDA0002125003080000082
Figure BDA0002125003080000083
Figure BDA0002125003080000084
Figure BDA0002125003080000085
其中，Pj Local(i)为图斑i对应等级j的空间局部概率，j＝1,2,3,…,5；
Wneighbor(i)＝max(Wij)
其中，图斑i局部最大影响权重对应的等级记为Rneighbor(i)，j＝1,2,3,…,5；
步骤4.2：将等级为j的图斑i的邻接图斑的数量和面积与全局中所有j等级图斑的邻接图斑数量和面积的最大值进行比较，得到数量比例PGlobal(Ai)与面积比例PGlobal(Bi)，计算图斑i的空间全局概率Pj Global(i)；
Figure BDA0002125003080000086
Figure BDA0002125003080000087
Figure BDA0002125003080000088
其中，j＝1,2,3,…,5。
步骤5：将相同研究区域内的植被覆盖度图及坡度图归一化后转换为栅格图；记x,y为栅格图行列号，z(x,y)为根据DEM数据得到的(x,y)像元处的高程值；结合水保措施情况，计算(x,y,z(x,y))像元的三维决策因子Pchange，Pchange越高更易发生水土流失；
Figure BDA0002125003080000091
Figure BDA0002125003080000092
Figure BDA0002125003080000093
其中，NDVI为植被指数，NDVImax、NDVImin为最大及最小NDVI值，V为植被覆盖度，L对应土地利用类型，Sx,y,z(x,y)为坡度，dz/dx和dz/dy分别表示像元及其邻域确定的x方向变化率和y方向变化率；W为水保措施情况，W＝0表示无水保措施，W＝1表示设置水保措施；
步骤6：根据“地理学第一定律”和空间关系，定义面向对象三维元胞的变化规则；
本实施例中，面向对象三维元胞的变化规则包括：
1、设定面向对象的选取规则并选择发生土壤侵蚀强度等级变化的图斑；
(4)如果图斑i的Wneighbor(i)＝1，该图斑的等级将变为Rneighbor(i)，记录其 IDcenter(i)于ChangeAll列表中，表示图斑内部全部发生等级变化；
(5)如果图斑i的Aj center(i)<a且Wneighbor(i)>b,该图斑内部等级全部变为 Rneighbor(i)，记录其IDcenter(i)于ChangeAll列表中；
(6)根据Pj Global(i)选取图斑内部部分区域发生等级变化的图斑，采用面积置信水平，数量置信水平及设置阈值三种方法进行计算；
将等级为j的所有图斑依据Pj Global(i)从大到小进行排序选取，面积选取为选取图斑面积和达到面积Aj choose时停止选取，数量选取为选取j等级图斑的设定比例数量的图斑，阈值选取为当Pj Global(i)大于设定阈值时进行选取；记录选取的图斑的IDcenter(i)于ChangePortion列表中表示图斑内部部分区域发生等级变化，公式如下：
Figure BDA0002125003080000094
其中，规则中a，b，c均代表置信水平，j＝1,2,3,…,5；
2、设定三维元胞自动机变化规则并改变土壤侵蚀强度等级；
将土壤侵蚀强度等级矢量图分别以图斑ID及图斑等级Rank为属性值分别转换为ID栅格及等级Rank栅格，将土地利用类型图转换为栅格；
遍历栅格，记录当前中心元胞为(x,y,z(x,y))，其中x,y为行列号，z(x,y)为根据DEM数据得到的(x,y)处高程值；t时刻中心元胞的ID为IDx,y,z(x,y)(t)，与其所位于的图斑i的ID号IDcenter(i)对应相等，土地利用类型为Lx,y,z(x,y)(t)，元胞的等级状态为Rx,y,z(x,y)(t)；以3×3为邻域定义x,y方向上的邻居元胞；其中，设定所有图斑 i的初始改变状态status(i)为0；
三维元胞自动机变化规则如下：
(1)当元胞(x,y,z(x,y))的IDx,y,z(x,y)(t)位于ChangeAll列表中，且三维决策因子Pchange不为0，即无水保措施，将该元胞的t+1时刻的等级Rx,y,z(x,y)(t+1)改变为 Rneighbor(i)；
(2)当元胞(x,y,z(x,y))与邻居元胞ID不相等时，即IDx,y,z(x,y)(t)不等于(x,y-1)位置的IDx,y-1,z(x,y-1)(t)或(x,y+1)位置的IDx,y+1,z(x,y+1)(t)时，元胞位于图斑的边界线；若IDx,y,z(x,y)(t)位于ChangePortion列表中且其邻域中存在邻居元胞ID号为IDneighbor(i) 时，记录像元位置为(m,n)，并以此为起始点像元开始进行局部变化；请见图7(A)，填充方式不同的区域代表不同的图斑，IDneighbor(i)为IDcenter(i)图斑的最大影响邻接图斑ID号；
局部等级变化规则有如下4种方法，请见图7(B)，表示按规则选取邻居元胞方向并在IDcenter(i)图斑的内部进行等级变化：
①在邻域内每次随机选择一个邻居元胞，依次生成连续成片的像元并将元胞 t+1时刻等级更改为Rneighbor(i)；
②选取邻居元胞中Pchange最大值更改等级为Rneighbor(i)；
③每次随机生成一个方向，如果其Pchange大于置信水平d，选中该元胞发生等级变化；若无邻居元胞符合条件，随机抽中一方向继续进行变化；
④计算邻居元胞Pchange占邻域范围所有邻居元胞Pchange和的百分比，即概率，以此为随机抽取的选中概率进行选取并改变等级；
选取过程中保证像元ID与起始点ID相等，即在该图斑内部发生变化，直到改变像元的面积和达到随机面积Achange(i)，设置status(i)为1；返回并由位置 (x,y+1)处像元开始继续遍历，公式如下：
Figure BDA0002125003080000111
其中，e为(0,1)之间的随机数,j＝1,2,3,…,5。
在步骤6的规则中，面向对象的选取规则为三维元胞自动机规则提供了决策依据，以达到综合考虑元胞局部空间信息、面向对象的区域空间信息和全局空间信息的目的，符合“地理学第一定律”。在大尺度研究下，高程值z(x,y)可视为几乎不发生改变。
步骤7：对研究区域进行水文分析；
本实施例中，根据DEM数据，确定元胞(x,y,z(x,y))位置处的流向Dirx,y,z(x,y)及流量Accx,y,z(x,y)；
Figure BDA0002125003080000112
Accx,y,z(x,y)＝F(Dirx,y,z(x,y),weight)
其中，n＝1,2,4,8,16,32,64,128，代表8个邻居元胞；dn表示元胞间距离，当 n＝1,4,16,64时，dn＝1；当n＝2,8,32,128时，dn＝1.5；流量函数F中Dirx,y,z(x,y)和 weight代表流向栅格和权重矩阵。
步骤8：基于发生等级变化后的土壤侵蚀强度等级图，重新计算空间概率，选取发生水土流失的区域。
本实施例中，选取发生水土流失的区域，选取规则为对于j等级的所有图斑，以局部最大影响等级与局部最大影响权重为双重条件进行降序排序，图斑i的 Rneighbor(i)与Wneighbor(i)越高，越易发生水土流失；按所排顺序选取图斑，直到选取的图斑面积和达到步骤3中计算的Aj happen为止；对步骤7得到的水文分析结果中提取出发生水土流失的区域，完成研究区域内水土流失模拟过程；其中， Aj happen是由贝叶斯原理对发生水土流失的面积进行分配得出各等级j发生水土流失区域的面积。
请见图2，图斑的ID号为i；面积Areai为图斑i的面积；Pi Local为图斑i的空间局部概率信息；Pi Global为图斑i的空间全局概率信息；Wi neighbor为图斑i的局部最大影响权重；ChangeAll、ChangePortion、NotChange分别记录全部发生、部分发生及不发生等级变化的图斑对应的ID号。
请见图3，在三维空间中，以研究的中心元胞为坐标原点，定义三维笛卡尔坐标系。每个元胞周围邻接26个邻居元胞。
请见图4，在数字高程模型(Digital Elevation Model，简称DEM)图中，以中心元胞定义图3所示的三维笛卡尔坐标系。
请见图5，Xi,j表示中心元胞，i、j为元胞对应像元在栅格图中的行列号。其在x,y方向上具有8个邻居元胞。
请见图6，n为流向代码。n＝1,2,4,8,16,32,64,128,代表8个邻居元胞。dn表示元胞间距离，dn＝1(n＝1,4,16,64),dn＝1.5(n＝2,8,32,128)。
请见图7，在栅格图中，填充方式不同的区域代表不同的图斑。图(A)图斑i 的ID号为IDcenter(i)，其最大影响邻接图斑ID号为IDneighbor(i)。若IDcenter(i)位于 ChangePortion列表中且其位于边界处的元胞(x,y)邻域中存在邻居元胞ID号为 IDneighbor(i)，进行图(B)的部分等级变化。以(x,y)为起始点，按规则选取邻居元胞方向并在该图斑内部进行等级变化。
应当理解的是，本说明书未详细阐述的部分均属于现有技术。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (7)

1.一种利用面向对象的三维元胞自动机进行水土流失模拟的方法，其特征在于，包括以下步骤：
步骤1：计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图；
步骤2：遍历土壤侵蚀强度等级矢量图中所有图斑，得到ID号为IDcenter(i)、等级为j的图斑i的面积Aj center(i)；记录等级为j的所有图斑的面积和为Aj rank；
Figure FDA0002125003070000011
其中，j＝1,2,3,…,N，N为等级总数；
步骤3：预测研究区域内发生水土流失的总面积Asum并进行分配得到各等级发生水土流失区域的面积，其中研究区域总面积为Aarea；
步骤4：遍历土壤侵蚀强度等级矢量图所有图斑，获取图斑空间联合概率；
步骤5：将相同研究区域内的植被覆盖度图及坡度图归一化后转换为栅格图；记x,y为栅格图行列号，z(x,y)为根据DEM数据得到的(x,y)像元处的高程值；结合水保措施情况，计算(x,y,z(x,y))像元的三维决策因子Pchange，Pchange越高更易发生水土流失；
Figure FDA0002125003070000012
Figure FDA0002125003070000013
Figure FDA0002125003070000014
其中，NDVI为植被指数，NDVImax、NDVImin为最大及最小NDVI值，V为植被覆盖度，L对应土地利用类型，Sx,y,z(x,y)为坡度，dz/dx和dz/dy分别表示像元及其邻域确定的x方向变化率和y方向变化率；W为水保措施情况，W＝0表示无水保措施，W＝1表示设置水保措施；
步骤6：根据“地理学第一定律”和空间关系，定义面向对象三维元胞的变化规则；
步骤7：对研究区域进行水文分析；
步骤8：基于发生等级变化后的土壤侵蚀强度等级图，重新计算空间概率，选取发生水土流失的区域。
2.根据权利要求1所述的利用面向对象的三维元胞自动机进行水土流失模拟的方法，其特征在于：步骤1中，利用三因子方法计算土壤侵蚀强度等级并得到土壤侵蚀强度等级矢量图，等级数目为5；其中，所述三因子为地形坡度、地表覆盖度和土地利用类型。
3.根据权利要求1所述的利用面向对象的三维元胞自动机进行水土流失模拟的方法，其特征在于：步骤3中，根据先验概率P(A)预测研究区域内发生水土流失的总面积Asum，其中研究区域总面积为Aarea；
由贝叶斯原理对发生水土流失的面积进行分配得出各等级j发生水土流失区域的面积Aj happen；
Figure FDA0002125003070000021
Asum＝Aarea×P(A)
Figure FDA0002125003070000022
Figure FDA0002125003070000023
其中，Aj rank为等级为j的所有图斑的面积和，j＝1,2,3,…,5；
Figure FDA0002125003070000024
为先验发生水土流失的面积比例，即先验概率；
P(Bj)表示发生第j等级水土流失的比例；
Figure FDA0002125003070000025
为先验发生水土流失时第j等级发生水土流失的面积比例。
4.根据权利要求1所述的利用面向对象的三维元胞自动机进行水土流失模拟的方法，其特征在于，步骤4的具体实现包括以下子步骤：
步骤4.1：计算图斑i对应的邻接j等级图斑占其所有邻接图斑的数量比例PLocal(Aij)与面积比例PLocal(Bij)，根据“地理学第一定律”，得到图斑i对应等级j的空间局部概率Pj Local(i)、局部影响权重Wij，得出图斑i的局部最大影响等级Rneighbor(i)及局部最大影响权重Wneighbor(i)；同时记录图斑i邻接图斑中等级为Rneighbor(i)且面积最大的图斑对应的ID号为IDneighbor(i)，表示图斑i的最大影响邻接图斑；
Figure FDA0002125003070000031
Figure FDA0002125003070000032
Figure FDA0002125003070000033
Figure FDA0002125003070000034
其中，Pj Local(i)为图斑i对应等级j的空间局部概率，j＝1,2,3,…,5；
Wneighbor(i)＝max(Wij)
其中，图斑i局部最大影响权重对应的等级记为Rneighbor(i)，j＝1,2,3,…,5；
步骤4.2：将等级为j的图斑i的邻接图斑的数量和面积与全局中所有j等级图斑的邻接图斑数量和面积的最大值进行比较，得到数量比例PGlobal(Ai)与面积比例PGlobal(Bi)，计算图斑i的空间全局概率Pj Global(i)；
Figure FDA0002125003070000035
Figure FDA0002125003070000036
Figure FDA0002125003070000037
其中，j＝1,2,3,…,5。
5.根据权利要求4所述的利用面向对象的三维元胞自动机进行水土流失模拟的方法，其特征在于，步骤6中所述面向对象三维元胞的变化规则，包括：
1、设定面向对象的选取规则并选择发生土壤侵蚀强度等级变化的图斑；
(1)如果图斑i的Wneighbor(i)＝1，该图斑的等级将变为Rneighbor(i)，记录其IDcenter(i)于ChangeAll列表中，表示图斑内部全部发生等级变化；
(2)如果图斑i的Aj center(i)<a且Wneighbor(i)>b,该图斑内部等级全部变为Rneighbor(i)，记录其IDcenter(i)于ChangeAll列表中；
(3)根据Pj Global(i)选取图斑内部部分区域发生等级变化的图斑，采用面积置信水平，数量置信水平及设置阈值三种方法进行计算；
将等级为j的所有图斑依据Pj Global(i)从大到小进行排序选取，面积选取为选取图斑面积和达到面积Aj choose时停止选取，数量选取为选取j等级图斑的设定比例数量的图斑，阈值选取为当Pj Global(i)大于设定阈值时进行选取；记录选取的图斑的IDcenter(i)于ChangePortion列表中表示图斑内部部分区域发生等级变化，公式如下：
Figure FDA0002125003070000041
其中，规则中a，b，c均代表置信水平，j＝1,2,3,…,5；
2、设定三维元胞自动机变化规则并改变土壤侵蚀强度等级；
将土壤侵蚀强度等级矢量图分别以图斑ID及图斑等级Rank为属性值分别转换为ID栅格及等级Rank栅格，将土地利用类型图转换为栅格；
遍历栅格，记录当前中心元胞为(x,y,z(x,y))，其中x,y为行列号，z(x,y)为根据DEM数据得到的(x,y)处高程值；t时刻中心元胞的ID为IDx,y,z(x,y)(t)，与其所位于的图斑i的ID号IDcenter(i)对应相等，土地利用类型为Lx,y,z(x,y)(t)，元胞的等级状态为Rx,y,z(x,y)(t)；以3×3为邻域定义x,y方向上的邻居元胞；其中，设定所有图斑i的初始改变状态status(i)为0；
三维元胞自动机变化规则如下：
(1)当元胞(x,y,z(x,y))的IDx,y,z(x,y)(t)位于ChangeAll列表中，且三维决策因子Pchange不为0，即无水保措施，将该元胞的t+1时刻的等级Rx,y,z(x,y)(t+1)改变为Rneighbor(i)；
(2)当元胞(x,y,z(x,y))与邻居元胞ID不相等时，即IDx,y,z(x,y)(t)不等于(x,y-1)位置的IDx,y-1,z(x,y-1)(t)或(x,y+1)位置的IDx,y+1,z(x,y+1)(t)时，元胞位于图斑的边界线；若IDx,y,z(x,y)(t)位于ChangePortion列表中且其邻域中存在邻居元胞ID号为IDneighbor(i)时，记录像元位置为(m,n)，并以此为起始点像元开始进行局部变化；
局部等级变化规则有如下4种方法：
①在邻域内每次随机选择一个邻居元胞，依次生成连续成片的像元并将元胞t+1时刻等级更改为Rneighbor(i)；
②选取邻居元胞中Pchange最大值更改等级为Rneighbor(i)；
③每次随机生成一个方向，如果其Pchange大于置信水平d，选中该元胞发生等级变化；若无邻居元胞符合条件，随机抽中一方向继续进行变化；
④计算邻居元胞Pchange占邻域范围所有邻居元胞Pchange和的百分比，即概率，以此为随机抽取的选中概率进行选取并改变等级；
选取过程中保证像元ID与起始点ID相等，即在该图斑内部发生变化，直到改变像元的面积和达到随机面积Achange(i)，设置status(i)为1；返回并由位置(x,y+1)处像元开始继续遍历，公式如下：
Figure FDA0002125003070000051
其中，e为(0,1)之间的随机数,j＝1,2,3,…,5。
6.根据权利要求5所述的利用面向对象的三维元胞自动机进行水土流失模拟的方法，其特征在于：步骤7中，根据DEM数据，确定元胞(x,y,z(x,y))位置处的流向Dirx,y,z(x,y)及流量Accx,y,z(x,y)；
Figure FDA0002125003070000052
Accx,y,z(x,y)＝F(Dirx,y,z(x,y),weight)
其中，n＝1,2,4,8,16,32,64,128，代表8个邻居元胞；dn表示元胞间距离，当n＝1,4,16,64时，dn＝1；当n＝2,8,32,128时，dn＝1.5；流量函数F中Dirx,y,z(x,y)和weight代表流向栅格和权重矩阵。
7.根据权利要求1-6任意一项所述的利用面向对象的三维元胞自动机进行水土流失模拟的***，其特征在于：步骤8中所述选取发生水土流失的区域，选取规则为对j等级的所有图斑，以局部最大影响等级与局部最大影响权重为双重条件进行降序排序，图斑i的Rneighbor(i)与Wneighbor(i)越高，越易发生水土流失；按所排顺序选取图斑，直到选取的图斑面积和达到步骤3中计算的Aj happen为止；对步骤7得到的水文分析结果中提取出发生水土流失的区域，完成研究区域内水土流失模拟过程；其中，Aj happen是由贝叶斯原理对发生水土流失的面积进行分配得出各等级j发生水土流失区域的面积。
CN201910619240.2A 2019-07-10 2019-07-10 利用面向对象的三维元胞自动机进行水土流失模拟的方法 Active CN110415346B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910619240.2A CN110415346B (zh) 2019-07-10 2019-07-10 利用面向对象的三维元胞自动机进行水土流失模拟的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910619240.2A CN110415346B (zh) 2019-07-10 2019-07-10 利用面向对象的三维元胞自动机进行水土流失模拟的方法

Publications (2)

Publication Number Publication Date
CN110415346A CN110415346A (zh) 2019-11-05
CN110415346B true CN110415346B (zh) 2022-11-25

Family

ID=68360882

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910619240.2A Active CN110415346B (zh) 2019-07-10 2019-07-10 利用面向对象的三维元胞自动机进行水土流失模拟的方法

Country Status (1)

Country Link
CN (1) CN110415346B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111047096A (zh) * 2019-12-16 2020-04-21 江西遥览科技有限公司 一种区域水土流失日预报算法
CN111858741A (zh) * 2020-07-22 2020-10-30 中国水利水电科学研究院 多期水土流失强度空间消长变化的可视化展示方法
CN113901640B (zh) * 2021-09-10 2024-05-28 中国矿业大学 一种基于元胞自动机的土壤养分流失仿真模拟方法
CN113505948B (zh) * 2021-09-13 2021-11-19 四川师范大学 基于贝叶斯网络的未来生态安全格局预测和优化方法
CN115358507B (zh) * 2022-06-30 2023-05-05 珠江水利委员会珠江流域水土保持监测中心站 生产建设项目扰动图斑水土流失风险识别评估方法
CN117049532B (zh) * 2023-10-11 2024-01-23 河北华运鸿业化工有限公司 一种制备固体氟化石墨的方法、***和设备

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE10325738B3 (de) * 2003-06-06 2004-12-02 Atlas Elektronik Gmbh Verfahren zur Generierung eines dreidimensionalen Geländemodells
CN103940974B (zh) * 2014-02-19 2016-03-09 西北农林科技大学 基于gis的中尺度流域土壤侵蚀时空动态分析方法
CN105004725B (zh) * 2015-08-04 2018-10-19 珠江水利委员会珠江水利科学研究院 一种水土保持综合治理土壤侵蚀变化量实时定量监测方法

Also Published As

Publication number Publication date
CN110415346A (zh) 2019-11-05

Similar Documents

Publication Publication Date Title
CN110415346B (zh) 利用面向对象的三维元胞自动机进行水土流失模拟的方法
CN105893972B (zh) 基于影像进行的违章建筑物自动监测方法及其实现***
Marjanović et al. Landslide susceptibility assessment using SVM machine learning algorithm
CN108830870B (zh) 基于多尺度结构学习的卫星影像高精度农田边界提取方法
CN110263111B (zh) 基于先验知识的土地利用/覆被信息时空监测方法
CN107067405B (zh) 基于尺度优选的遥感影像分割方法
CN107977992A (zh) 一种基于无人机激光雷达的建筑物变化检测方法及装置
Saxena et al. Capturing heterogeneous urban growth using SLEUTH model
CN113192086B (zh) 地质灾害隐患变形强度分布图的生成方法和存储介质
Chang et al. Applying a modified VIKOR method to classify land subdivisions according to watershed vulnerability
CN111028255A (zh) 基于先验信息与深度学习的农田区域预筛选方法及装置
CN107330422A (zh) 一种基于高精度数字高程模型对半干旱地区进行微地形分类的方法
Wang et al. The isotropic organization of DEM structure and extraction of valley lines using hexagonal grid
CN102073867B (zh) 一种遥感图像分类方法及装置
CN109658477A (zh) 一种基于lidar数据的dem生成算法
CN111611960B (zh) 一种基于多层感知神经网络大区域地表覆盖分类方法
US11544297B2 (en) Computer-based method and system for urban planning
CN110059860A (zh) 一种城市的公共充电站位置布局多目标优化方法
CN106780586A (zh) 一种基于地面激光点云的太阳能潜力评估方法
Xu et al. Building height calculation for an urban area based on street view images and deep learning
Ahmad et al. The cellular automata approach in dynamic modelling of land use change detection and future simulations based on remote sensing data in Lahore Pakistan
CN114723121A (zh) 基于gis的野外复杂地形路径规划方法
Vaidya et al. Classifying heterogeneous urban form into local climate zones using supervised learning and greedy clustering incorporating Landsat dataset
CN108921194B (zh) 一种地形坡位分类自适应的聚类方法
CN113705326B (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