CN112950779B - 一种度量地貌破碎程度的栅格化曲面的构建方法及系统 - Google Patents
一种度量地貌破碎程度的栅格化曲面的构建方法及系统 Download PDFInfo
- Publication number
- CN112950779B CN112950779B CN202110221023.5A CN202110221023A CN112950779B CN 112950779 B CN112950779 B CN 112950779B CN 202110221023 A CN202110221023 A CN 202110221023A CN 112950779 B CN112950779 B CN 112950779B
- Authority
- CN
- China
- Prior art keywords
- landform
- curved surface
- rasterized
- unit
- data
- 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
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
- G06T7/00—Image analysis
- G06T7/60—Analysis of geometric attributes
- G06T7/64—Analysis of geometric attributes of convexity or concavity
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Geometry (AREA)
- Software Systems (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Remote Sensing (AREA)
- Computer Graphics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Processing Or Creating Images (AREA)
- Instructional Devices (AREA)
Abstract
本发明涉及栅格化地形曲面构建技术领域，涉及一种度量地貌破碎程度的栅格化曲面的构建方法及系统，其包括以下步骤：一、计算机获取地貌模型；二、进行预处理，得到预处理后数据；三、进行窗口标准差统计，得到表征地貌高程差异程度的栅格化曲面数据；四、进行地貌高差提取，得到表征地貌隆起与水平切割深度的栅格化曲面数据；五、进行沟壑密度提取，得到表征地形切割密度的栅格化曲面数据；六、计算得到度量地貌破碎程度的栅格化地形曲面。本发明能够准确描绘和表达地貌破碎度的真实程度，有利于构建各种地学指标，刻画相应的地理环境形态。
Description
技术领域
本发明涉及栅格化地形曲面构建技术领域，进一步的说，涉及用栅格化电子地图模拟表达、虚拟和计算地球资源和环境的时空分布的技术领域，具体地说，涉及一种度量地貌破碎程度的栅格化曲面的构建方法及系统。
背景技术
自然地理环境中，一定空间尺度上的土地单元，地貌形态并不是完整、有规律的分布，而是具有一定数量的破碎小斑块，呈连续或不连续分布特征。地貌破碎与连续不断升高的巨大山体地表风的方向具有一致性，破碎山地的风向改变更为频繁，能导致特殊气流变化，影响温度分布、水汽输送、局地风场及风蚀、水蚀等地表过程，进而影响土壤分布及成分形成。地貌破碎也是水平地带性与垂直地带性等地域分异规律形成的主要驱动因子，其水平方向上揭示研究区从海岸沿线到山体水平带谱的分布状况，垂直方向上揭示从山麓到山顶的地貌垂直变化特征。构建度量地貌破碎的栅格化地形曲面对于山体分布密度大且变化复杂的山地区域具有显著的科学及实践意义，诸如我国南部云贵高原、云岭地区。地貌破碎程度取决于一定地理范围内同级别地貌单元内呈破碎状分布的小斑块数量的多少。数字高程模型DEM(Digital Elevation Model)和数字地形模型DTM(Digital Terrain Model)用栅格化高程(海拔)数据模拟表达地表的高低起伏或者陡峭。在某些地理指标的应用场合，直接使用DEM或DTM，难以直接描述地貌的破碎化程度。
发明内容
本发明的内容是提供一种度量地貌破碎程度的栅格化曲面的构建方法及系统，其能够克服现有技术的某种或某些缺陷。
根据本发明的一种度量地貌破碎程度的栅格化曲面的构建方法，其包括以下步骤：
一、计算机获取地貌模型D；
二、对所述地貌模型D进行预处理，得到预处理后数据D’；
三、对所述预处理后数据D’进行窗口标准差统计，得到整个地理区域中不同地貌单元内的、在水平方向上的高程差异程度的栅格化曲面数据，记为Std；
四、对所述预处理后数据进行地貌高差提取，得到整个地理区域中不同地貌单元内的、在垂直方向上表征的地貌隆起与水平切割深度的栅格化曲面数据，记为Terr；
五、对所述预处理后数据进行沟壑密度提取，得到整个地理区域中不同地貌单元内表征地形切割密度的栅格化曲面数据，记为M；
六、使用公式A＝Std×Terr×M，得到度量地貌破碎程度的栅格化地形曲面，记为A。
作为优选，步骤一中，所述地貌模型D采用计算机表达的数字高程模型DEM或者数字地形模型DTM。
作为优选，步骤二中，预处理的方法为：对所述地貌模型数据D，进行格式转换、投影变换和噪声滤波处理，得到预处理后数据D’；所述投影变换将地理坐标投影变换为以米为长度单位的平面坐标和高程坐标。
作为优选，步骤三中，窗口标准差统计的方法为：将所述预处理后数据D’，以3*3大小的栅格分析窗口为单元，使用空间分析焦点统计功能，统计各单元的标准差值，得到整个地理区域中不同地貌单元标准差值，生成表征地貌高程差异程度的栅格化曲面数据Std。
作为优选，步骤四中，地貌高差提取的方法为：
a、将所述预处理后的数据D’，进行分水岭提取，得到覆盖整个地理区域各分水岭地貌单元多边形F_vi，其中i为各分水岭地貌单元的标识号，i＝1……m，m为所述各分水岭地貌单元的总个数；
b、确定各所述多边形F_vi的最大内接椭圆Yi和各所述最大内接椭圆的圆心坐标ki、li；
c、根据各所述圆心坐标，使用公式ui(ki,li)＝(h1+h2+...+hj)/m，确定各所述多边形内的地貌模型栅格点的高程算术平均值，其中ui(ki,li)为圆心位置的数字高程的算术平均值求算，hj为每个栅格点的数字高程值，m为多边形内全部栅格点的数量；
d、根据各所述圆心坐标和所述算术平均值采用克里金插值法，得到地貌趋势面W；
e、根据所述地貌模型D’和所述地貌趋势面W，将所述地貌模型D’减去所述地貌趋势面W，得到地貌相对高度度量Terr_o；
f、将所述地貌相对高度度量Terr_o使用离差标准化法进行归一化处理，得到度量地貌隆起与水平切割深度的栅格化曲面数据Terr。
作为优选，步骤五中，沟壑密度提取的方法为：
1)、调入所述预处理后的数据D’，将每个单元海拔值与周边8个单元比较，确定当前单元水流的最大下降方向；
2)、累加下降到该点的水文汇流量，储存为F_accum；
3)、人机交互确定阈值a，F_accum与a进行大小比较，大于a的单元判断为沟谷地貌，单元值记为1，否则记为Nodata，储存为E，对E矢量化得到整个地理区域上沟谷线数据E_v；
4)、根据各所述分水岭地貌单元多边形F_vi与沟谷线数据E_v，进行空间标识，得到表征各分水岭内河流分布的沟谷线数据E_vi,j，其中j为整个地理区域上第i个分水岭地貌单元多边形F_vi内各沟谷线的标识号，j＝1……n，n为所述各分水岭地貌单元内沟谷线的总条数；
5)、根据所述F_vi及E_vi,j数据，进行空间叠加分析和几何特征运算，得到各分水岭地貌单元内沟谷线的总长度Li与所对应的分水岭面积Si；
6)、将所述各分水岭地貌单元内沟谷线的总长度Li除以所述所对应的分水岭面积Si，得到沟壑密度值Mi，公式为Mi＝Li/Si；
7)、根据所述整个地理区域上各分水岭地貌单元多边形F_vi和所述沟壑密度值Mi，进行矢量转栅格处理，得到整个地理区域中表征地形切割密度的栅格化曲面数据M。
作为优选，步骤六中，度量地貌破碎程度的栅格化地形曲面A的的获得方法为：将所述表征地貌高程变化差异程度的栅格化曲面数据Std、所述表征地貌隆起与水平切割深度的栅格化曲面数据Terr和所述表征地形切割密度的栅格化曲面数据M，使用公式A＝Std×Terr×M，进行栅格乘积运算，得到度量地貌破碎程度的栅格化地形曲面A。
本发明还提供了一种度量地貌破碎程度的栅格化曲面的构建系统，其采用上述的一种度量地貌破碎程度的栅格化曲面的构建方法，并包括：
地貌数据模块，用于获取地貌模型；
预处理模块，用于对所述地貌模型进行预处理，得到预处理后数据；
地貌变化差异曲面计算模块，用于对所述预处理后数据进行窗口标准差统计，得到整个地理区域中不同地貌单元内的、在水平方向上的高程差异程度的栅格化曲面数据；
地貌隆起与水平切割深度曲面计算模块，用于对所述预处理后数据进行地貌高差提取，得到整个地理区域中不同地貌单元内的、在垂直方向上表征的地貌隆起与水平切割深度的栅格化曲面数据；
地形切割密度曲面计算模块，用于对所述预处理后数据进行沟壑密度提取，得到整个地理区域中不同地貌单元内表征地形切割密度的栅格化曲面数据；
地貌破碎程度曲面计算模块，用于根据所述地貌高程变化差异程度的栅格化曲面数据、所述地貌隆起与水平切割深度的栅格化曲面数据、所述地形切割密度的栅格化曲面数据，得到度量地貌破碎程度的栅格化地形曲面。
作为优选，预处理模块包括预处理单元，预处理单元用于对所述地貌数据进行格式转换、投影变换和噪声滤波处理，得到预处理后数据；所述投影变换将地理坐标投影变换为以米为长度单位的平面坐标；
地貌变化差异曲面计算模块包括地貌变化差异曲面确定单元，地貌变化差异曲面确定单元用于将所述预处理后数据D’，以3*3大小的栅格分析窗口为单元，使用空间分析焦点统计功能，统计各单元的标准差值，得到整个地理区域中不同地貌单元标准差值，生成表征地貌高程变化差异程度的栅格化曲面数据Std；
地貌隆起与水平切割深度曲面计算模块包括分水岭提取单元、内接椭圆确定单元、高程算术平均值确定单元、地貌趋势面确定单元、地貌相对高度曲面确定单元和地貌隆起与水平切割深度曲面确定单元；
分水岭提取单元用于对所述预处理后数据D’进行分水岭提取，得到各分水岭地貌单元多边形F_vi，其中i为各分水岭地貌单元的标识号，i＝1……m，m为所述各分水岭地貌单元的总个数；
内接椭圆确定单元用于确定各分水岭地貌单元多边形F_vi的最大内接椭圆Yi和各所述最大内接椭圆的圆心坐标ki、li；
高程算术平均值确定单元用于根据各所述圆心坐标ki、li，确定各所述多边形内的地貌模型栅格点的高程算术平均值ui(ki,li)；
地貌趋势面确定单元用于根据各所述圆心坐标ki、li和所述算术平均值ui(ki,li)采用克里金插值法，得到地貌趋势面W；
地貌相对高度曲面确定单元用于根据所述地貌模型D’和所述地貌趋势面W，将所述地貌模型D’减去所述地貌趋势面W，得到地貌相对高度度量Terr_o；
地貌隆起与水平切割深度曲面确定单元用于将所述地貌相对高度度量Terr_o使用离差标准化法进行归一化处理，得到度量地貌隆起与水平切割深度的栅格化曲面数据Terr；
地形切割密度曲面计算模块包括沟谷线提取单元和地形切割密度曲面确定单元；
沟谷线提取单元用于调入所述预处理后的数据D’，将每个单元海拔值与周边8个单元比较，确定当前单元水流的最大下降方向；累加下降到该点的水文汇流量，储存为F_accum；人机交互确定阈值a，F_accum与a进行大小比较，大于a的单元判断为沟谷地貌，单元值记为1，否则记为Nodata，储存为E，对E矢量化得到整个地理区域上沟谷线数据E_v；根据各所述分水岭地貌单元多边形F_vi与沟谷线数据E_v，进行空间标识，得到表征各分水岭内河流分布的沟谷线数据E_vi,j，其中j为整个地理区域上第i个分水岭地貌单元多边形F_vi内各沟谷线的标识号，j＝1……n，n为所述各分水岭地貌单元内沟谷线的总条数；
地形切割密度曲面确定单元用于根据所述F_vi及E_vi,j数据，进行空间叠加分析和几何特征运算，得到各分水岭地貌单元内沟谷线的总长度Li与所对应的分水岭面积Si；将所述各分水岭地貌单元内沟谷线的总长度Lj除以所述所对应的分水岭面积Si，得到沟壑密度值Mi，公式为Mi＝Li/Si；根据所述整个地理区域上各分水岭地貌单元多边形边界数据F_v和所述沟壑密度值Mi，进行矢量转栅格处理，得到整个地理区域中表征地形切割密度的栅格化曲面数据M。
作为优选，地貌破碎程度曲面计算模块包括地貌破碎程度曲面确定单元，地貌破碎程度曲面确定单元用于将所述地貌高程差异程度的栅格化曲面数据Std、所述地貌隆起与水平切割深度的栅格化曲面数据Terr和所述地形切割密度的栅格化曲面数据M，使用公式A＝Std×Terr×M，进行栅格乘积运算，得到度量地貌破碎程度的栅格化地形曲面A。
本发明能够准确描绘和表达地貌破碎度的真实程度，有利于构建各种地学指标，刻画相应的地理环境形态。本发明用以表征地貌多样性，用栅格大数据模拟表达和计算地球资源变量的多样性；并构建相应的计算机数字化处理系统。
附图说明
图1为实施例1中一种度量地貌破碎度的栅格化地形曲面的构建方法的流程图；
图2为实施例2中一种度量地貌破碎度的栅格化地形曲面的构建系统的结构框图。
具体实施方式
为进一步了解本发明的内容，结合附图和实施例对本发明作详细描述。应当理解的是，实施例仅仅是对本发明进行解释而并非限定。
实施例1
如图1所示，本实施例提供了一种度量地貌破碎程度的栅格化曲面的构建方法，其包括以下步骤：
一、计算机获取地貌模型D；利用数字地形模型，进行栅格重采样或地形尺度变换，得到一定像元尺度的地貌表达数据，得到所述地貌模型D。所述地貌模型D采用计算机表达的数字高程模型DEM或者数字地形模型DTM。
二、对所述地貌模型D进行预处理，得到预处理后数据D’；具体为：对所述地貌模型数据D，进行格式转换、投影变换和噪声滤波处理，得到预处理后数据D’；所述投影变换使用地理信息系统软件的投影算法将具有地理坐标的地貌模型数据投影变换为以米为长度单位的平面坐标和高程坐标。所述投影算法为一般的地理信息系统常规算法。
三、对所述预处理后数据D’进行窗口标准差统计计算，得到整个地理区域中不同地貌单元内的、在水平方向(X-Y平面)上的高程差异程度的栅格化曲面数据，记为Std；具体为：利用地理信息系统软件栅格像元统计算法，以所述预处理后数据D’为分析对象，以3*3像元大小的矩形分析窗口为统计单元，使用窗口滑动法遍历整个地理区域各像元，统计各单元的标准差值，得到整个地理区域中不同单元的变化差异，生成表征地貌高程变化差异程度的栅格化曲面数据Std。用于矩形窗口像元统计的算法是一般地理信息系统公开的常规算法，称为焦点统计。
四、对所述预处理后数据进行地貌高差提取，得到整个地理区域中不同地貌单元内的、在垂直方向上表征的地貌隆起与水平切割深度的栅格化曲面数据，记为Terr；具体为：
a、以所述预处理后的数据D’为计算机输入数据，使用地理信息系统软件中的分水岭提取算法，得到覆盖整个地理区域各分水岭地貌单元多边形F_vi，其中i为各分水岭地貌单元的标识号，i＝1……m，m为所述各分水岭地貌单元的总个数；
b、使用地理信息系统软件中的最大内接椭圆和圆心坐标的计算算法，确定各所述多边形F_vi的最大内接椭圆Yi和各所述最大内接椭圆的圆心坐标ki、li；以圆心坐标的坐标位置作为该流域即该多边形的地貌趋势面的代表位置；
c、根据各所述圆心坐标，使用公式ui(ki,li)＝(h1+h2+...+hj)/m，确定各所述多边形内的地貌模型栅格点的高程算术平均值，其中ui(ki,li)为圆心位置的数字高程的算术平均值求算，hj为每个栅格点的数字高程值，m为多边形内全部栅格点的数量；
d、根据各所述圆心坐标和所述算术平均值采用克里金插值法，得到地貌趋势面W；
e、根据所述地貌模型D’和所述地貌趋势面W，将所述地貌模型D’减去所述地貌趋势面W，得到地貌相对高度度量Terr_o；
f、将所述地貌相对高度度量Terr_o使用离差标准化法进行归一化处理，得到度量地貌隆起与水平切割深度的栅格化曲面数据Terr。
分水岭提取、最大内接椭圆和圆心坐标的计算算法，是一般的地理信系统的公开的常规算法。离差标准化法是一般的数理统计公开的常规方法。
五、对所述预处理后数据进行沟壑密度提取，得到整个地理区域中不同地貌单元内表征地形水平切割密度的栅格化曲面数据，记为M；具体为：
1)、以所述预处理后的数据D’为计算机输入数据，使用地理信息系统软件水文分析中流向计算算法，将每个单元海拔值与周边8个单元比较，确定当前单元水流的最大下降方向；
2)、使用地理信息系统软件水文分析中汇流量计算算法，计算累加下降到该点的水文汇流量，储存为F_accum；
3)、人机交互确定阈值a，使用地理信息系统软件条件选取及重分类算法，将F_accum与a进行大小比较，大于a的单元判断为沟谷地貌，单元值记为1，否则记为Nodata，储存为E，使用地理信息系统软件栅格数据转矢量数据算法，对E矢量化得到整个地理区域上沟谷线数据E_v；
4)、使用地理信息系统软件空间连接算法，根据各所述分水岭地貌单元多边形F_vi与沟谷线数据E_v，进行空间标识，得到表征各分水岭内河流分布的沟谷线数据E_vi,j，其中j为整个地理区域上第i个分水岭地貌单元多边形F_vi内各沟谷线的标识号，j＝1……n，n为所述各分水岭地貌单元内沟谷线的总条数；
5)、根据所述F_vi及E_vi,j数据，进行空间叠加分析和几何特征运算，得到各分水岭地貌单元内沟谷线的总长度Li与所对应的分水岭面积Si；
6)、将所述各分水岭地貌单元内沟谷线的总长度Li除以所述所对应的分水岭面积Si，得到沟壑密度值Mi，公式为Mi＝Li/Si；
7)、根据所述整个地理区域上各分水岭地貌单元多边形F_vi和所述沟壑密度值Mi，进行矢量转栅格处理，得到整个地理区域中表征地形切割密度的栅格化曲面数据M。
上述水文分析中流向计算算法、汇流量计算算法、条件选取及重分类算法、栅格数据转矢量数据算法、空间连接算法是一般的地理信息系统公开的常规算法。
六、使用公式A＝Std×Terr×M，得到度量地貌破碎程度的栅格化地形曲面，记为A；具体为：将所述地貌高程差异程度的栅格化曲面数据Std、所述地貌隆起与水平切割深度的栅格化曲面数据Terr和所述地形切割密度的栅格化曲面数据M为计算机输入数据，使用地理信息系统软件栅格运算算法，使用公式A＝Std×Terr×M，进行栅格乘积运算，得到度量地貌破碎程度的栅格化地形曲面A。
栅格运算算法是一般的地理信系统的公开的常规算法。
如果空间上一片连接区域的地貌破碎度小于3，表示该地理区域为小地貌分布密度较小的区域；如果在[3,9]范围间，表示该地理区域为小地貌分布密度中等的区域；如果大于9，表示该地理区域为小地貌分布密度较大的区域。
利用上述的地貌破碎度，或地貌破碎度和不同大小的栅格单元面积等，即可得到某山地区域单位平面面积上地形地貌的破碎程度，进而可以用于复杂山区精细化气候模拟及其他地学研究。
以上方法涉及的DEM或DTM、D’、Std、Terr、M、A，是计算机中用地图坐标系统存储、管理和运算的数字矩阵，按数据格式称为栅格数据；矩阵的行列坐标、地图坐标可以进行相互转换。以上方法的实施，可以借助商业化地理信系统软件工具，在计算机中运算和实现。
实施例2
如图2所示，本实施例提供了一种度量地貌破碎程度的栅格化曲面的构建系统，其采用实施例1中的一种度量地貌破碎程度的栅格化曲面的构建方法，并包括：
地貌数据模块，用于获取地貌模型；
预处理模块，用于对所述地貌模型进行预处理，得到预处理后数据；
地貌变化差异曲面计算模块，用于对所述预处理后数据进行窗口标准差统计，得到整个地理区域中不同地貌单元内的、在水平方向(X-Y平面)上的高程差异程度的栅格化曲面数据；
地貌隆起与水平切割深度曲面计算模块，用于对所述预处理后数据进行地貌高差提取，得到整个地理区域中不同地貌单元内的、在垂直方向上表征的地貌隆起与水平切割深度的栅格化曲面数据；
地形切割密度曲面计算模块，用于对所述预处理后数据进行沟壑密度提取，得到整个地理区域中不同地貌单元内表征地形水平切割密度的栅格化曲面数据；
地貌破碎程度曲面计算模块，用于根据所述地貌高程差异程度的栅格化曲面数据、所述地貌隆起与水平切割深度的栅格化曲面数据、所述地形切割密度的栅格化曲面数据，得到度量地貌破碎程度的栅格化地形曲面。
地貌数据模块、预处理模块、地貌变化差异曲面计算模块、地貌隆起与水平切割深度曲面计算模块、地形切割密度曲面计算模块和地貌破碎程度曲面计算模块依次连接。
预处理模块包括预处理单元，预处理单元用于对所述地貌数据进行格式转换、投影变换和噪声滤波处理，得到预处理后数据；所述投影变换将地理坐标投影变换为以米为长度单位的平面坐标；
地貌变化差异曲面计算模块包括地貌变化差异曲面确定单元，地貌变化差异曲面确定单元用于将所述预处理后数据D’，以3*3大小的栅格分析窗口为单元，使用空间分析焦点统计功能，统计各单元的标准差值，得到整个地理区域中不同地貌单元标准差值，生成表征地貌高程变化差异程度的栅格化曲面数据Std；
地貌隆起与水平切割深度曲面计算模块包括分水岭提取单元、内接椭圆确定单元、高程算术平均值确定单元、地貌趋势面确定单元、地貌相对高度曲面确定单元和地貌隆起与水平切割深度曲面确定单元；
分水岭提取单元用于对所述预处理后数据D’进行分水岭提取，得到各分水岭地貌单元多边形F_vi，其中i为各分水岭地貌单元的标识号，i＝1……m，m为所述各分水岭地貌单元的总个数；
内接椭圆确定单元用于确定各分水岭地貌单元多边形F_vi的最大内接椭圆Yi和各所述最大内接椭圆的圆心坐标ki、li；
高程算术平均值确定单元用于根据各所述圆心坐标ki、li，确定各所述多边形内的地貌模型栅格点的高程算术平均值ui(ki,li)；
地貌趋势面确定单元用于根据各所述圆心坐标ki、li和所述算术平均值ui(ki,li)采用克里金插值法，得到地貌趋势面W；
地貌相对高度曲面确定单元用于根据所述地貌模型D’和所述地貌趋势面W，将所述地貌模型D’减去所述地貌趋势面W，得到地貌相对高度度量Terr_o；
地貌隆起与水平切割深度曲面确定单元用于将所述地貌相对高度度量Terr_o使用离差标准化法进行归一化处理，得到度量地貌隆起与水平切割深度的栅格化曲面数据Terr；
地形切割密度曲面计算模块包括沟谷线提取单元和地形切割密度曲面确定单元；
沟谷线提取单元用于调入所述预处理后的数据D’，将每个单元海拔值与周边8个单元比较，确定当前单元水流的最大下降方向；累加下降到该点的水文汇流量，储存为F_accum；人机交互确定阈值a，F_accum与a进行大小比较，大于a的单元判断为沟谷地貌，单元值记为1，否则记为Nodata，储存为E，对E矢量化得到整个地理区域上沟谷线数据E_v；根据各所述分水岭地貌单元多边形F_vi与沟谷线数据E_v，进行空间标识，得到表征各分水岭内河流分布的沟谷线数据E_vi,j，其中j为整个地理区域上第i个分水岭地貌单元多边形F_vi内各沟谷线的标识号，j＝1……n，n为所述各分水岭地貌单元内沟谷线的总条数；
地形切割密度曲面确定单元用于根据所述F_vi及E_vi,j数据，进行空间叠加分析和几何特征运算，得到各分水岭地貌单元内沟谷线的总长度Li与所对应的分水岭面积Si；将所述各分水岭地貌单元内沟谷线的总长度Lj除以所述所对应的分水岭面积Si，得到沟壑密度值Mi，公式为Mi＝Li/Si；根据所述整个地理区域上各分水岭地貌单元多边形边界数据F_v和所述沟壑密度值Mi，进行矢量转栅格处理，得到整个地理区域中表征地形切割密度的栅格化曲面数据M。
地貌破碎程度曲面计算模块包括地貌破碎程度曲面确定单元，地貌破碎程度曲面确定单元用于将所述地貌高程差异程度的栅格化曲面数据Std、所述地貌隆起与水平切割深度的栅格化曲面数据Terr和所述地形切割密度的栅格化曲面数据M，使用公式A＝Std×Terr×M，进行栅格乘积运算，得到度量地貌破碎程度的栅格化地形曲面A。
本实施例以DEM或DTM为基本数据，将数据调入计算机中，可显示和表达为地图坐标系统中计算机栅格图像。首先，计算机沿着数据集表达的高程和平面坐标，根据中心像元及邻域3*3窗口内所有像元具有标准差，计算表征地貌高程变化差异程度的栅格化曲面数据。其次，计算机沿着数据集表达的高程和平面坐标，经分水岭提取后，以每个分水岭多边形为单元，生成地势趋势面曲面；以DEM或DTM的像元为单位，计算表征地势高差的地貌隆起和切割深度曲面数据。再次，计算机沿着数据集表达的高程和平面坐标，经分水岭和沟谷线提取后，以每个分水岭多边形为单元，计算各分水岭多边形沟壑密度，从而获得表征地貌切割密度的曲面。最后，计算机按照地貌变化差异曲面、地貌隆起与水平切割深度曲面、地貌切割密度曲面的最小单元，沿着各数据集表达的平面坐标，逐像元做乘积运算，得到整个地理区域上度量地貌破碎程度的栅格化地形曲面。
地貌破碎度曲面是地图坐标系中的表达地貌破碎程度的栅格图像，可用于任何的地学分析中。其数值小于3的栅格点代表的为小地貌分布密度较小的区域；数字在3至9之间的栅格点代表小地貌分布密度中等的区域；数值大于9的栅格点小地貌分布密度较大的区域。
以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均应属于本发明的保护范围。
Claims (7)
1.一种度量地貌破碎程度的栅格化曲面的构建方法，其特征在于：包括以下步骤：
一、计算机获取地貌模型D；
二、对所述地貌模型D进行预处理，得到预处理后数据D’；
三、对所述预处理后数据D’进行窗口标准差统计，得到整个地理区域中不同地貌单元内的、在水平方向上的地貌高程差异程度的栅格化曲面数据，记为Std；
四、对所述预处理后数据进行地貌高差提取，得到整个地理区域中不同地貌单元内的、在垂直方向上表征的地貌隆起与水平切割深度的栅格化曲面数据，记为Terr；
地貌高差提取的方法为：
a、将所述预处理后的数据D’，进行分水岭提取，得到覆盖整个地理区域各分水岭地貌单元多边形F_vi，其中i为各分水岭地貌单元的标识号，i＝1……m，m为所述各分水岭地貌单元的总个数；
b、确定各所述多边形F_vi的最大内接椭圆Yi和各所述最大内接椭圆的圆心坐标ki、li；
c、根据各所述圆心坐标，使用公式ui(ki,li)＝(h1+h2+...+hj)/m，确定各所述多边形内的地貌模型栅格点的高程算术平均值，其中ui(ki,li)为圆心位置的数字高程的算术平均值求算，hj为每个栅格点的数字高程值，m为多边形内全部栅格点的数量；
d、根据各所述圆心坐标和所述算术平均值采用克里金插值法，得到地貌趋势面W；
e、根据所述地貌模型D’和所述地貌趋势面W，将所述地貌模型D’减去所述地貌趋势面W，得到地貌相对高度度量Terr_o；
f、将所述地貌相对高度度量Terr_o使用离差标准化法进行归一化处理，得到度量地貌隆起与水平切割深度的栅格化曲面数据Terr；
五、对所述预处理后数据进行沟壑密度提取，得到整个地理区域中不同地貌单元内表征地形切割密度的栅格化曲面数据，记为M；
沟壑密度提取的方法为：
1)、调入所述预处理后的数据D’，将每个单元海拔值与周边8个单元比较，确定当前单元水流的最大下降方向；
2)、累加下降到该点的水文汇流量，储存为F_accum；
3)、人机交互确定阈值a，F_accum与a进行大小比较，大于a的单元判断为沟谷地貌，单元值记为1，否则记为Nodata，储存为E，对E矢量化得到整个地理区域上沟谷线数据E_v；
4)、根据各所述分水岭地貌单元多边形F_vi与沟谷线数据E_v，进行空间标识，得到表征各分水岭内河流分布的沟谷线数据E_vi,j，其中j为整个地理区域上第i个分水岭地貌单元多边形F_vi内各沟谷线的标识号，j＝1……n，n为所述各分水岭地貌单元内沟谷线的总条数；
5)、根据所述F_vi及E_vi,j数据，进行空间叠加分析和几何特征运算，得到各分水岭地貌单元内沟谷线的总长度Li与所对应的分水岭面积Si；
6)、将所述各分水岭地貌单元内沟谷线的总长度Li除以所述所对应的分水岭面积Si，得到沟壑密度值Mi，公式为Mi＝Li/Si；
7)、根据所述整个地理区域上各分水岭地貌单元多边形F_vi和所述沟壑密度值Mi，进行矢量转栅格处理，得到整个地理区域中表征地形切割密度的栅格化曲面数据M；
六、使用公式A＝Std×Terr×M，得到度量地貌破碎程度的栅格化地形曲面，记为A。
2.根据权利要求1所述的一种度量地貌破碎程度的栅格化曲面的构建方法，其特征在于：步骤一中，所述地貌模型D采用计算机表达的数字高程模型DEM或者数字地形模型DTM。
3.根据权利要求2所述的一种度量地貌破碎程度的栅格化曲面的构建方法，其特征在于：步骤二中，预处理的方法为：对所述地貌模型数据D，进行格式转换、投影变换和噪声滤波处理，得到预处理后数据D’；所述投影变换将地理坐标投影变换为以米为长度单位的平面坐标和高程坐标。
4.根据权利要求3所述的一种度量地貌破碎程度的栅格化曲面的构建方法，其特征在于：步骤三中，窗口标准差统计的方法为：将所述预处理后数据D’，以3*3大小的栅格分析窗口为单元，使用空间分析焦点统计功能，统计各单元的标准差值，得到整个地理区域中不同地貌单元标准差值，生成表征地貌高程差异程度的栅格化曲面数据Std。
5.根据权利要求4所述的一种度量地貌破碎程度的栅格化曲面的构建方法，其特征在于：步骤六中，度量地貌破碎程度的栅格化地形曲面A的获得方法为：将所述地貌高程差异程度的栅格化曲面数据Std、所述地貌隆起与水平切割深度的栅格化曲面数据Terr和所述地形切割密度的栅格化曲面数据M，使用公式A＝Std×Terr×M，进行栅格乘积运算，得到度量地貌破碎程度的栅格化地形曲面A。
6.一种度量地貌破碎程度的栅格化曲面的构建系统，其特征在于：其采用如权利要求1-5中所述的任意一种度量地貌破碎程度的栅格化曲面的构建方法，并包括：
地貌数据模块，用于获取地貌模型；
预处理模块，用于对所述地貌模型进行预处理，得到预处理后数据；
地貌变化差异曲面计算模块，用于对所述预处理后数据进行窗口标准差统计，得到整个地理区域中不同地貌单元内的、在水平方向上的高程差异程度的栅格化曲面数据；
地貌隆起与水平切割深度曲面计算模块，用于对所述预处理后数据进行地貌高差提取，得到整个地理区域中不同地貌单元内的、在垂直方向上表征的地貌隆起与水平切割深度的栅格化曲面数据；
地形切割密度曲面计算模块，用于对所述预处理后数据进行沟壑密度提取，得到整个地理区域中不同地貌单元内表征地形切割密度的栅格化曲面数据；
地貌破碎程度曲面计算模块，用于根据所述地貌高程变化差异程度的栅格化曲面数据、所述地貌隆起与水平切割深度的栅格化曲面数据、所述地形切割密度的栅格化曲面数据，得到度量地貌破碎程度的栅格化地形曲面；
预处理模块包括预处理单元，预处理单元用于对所述地貌数据进行格式转换、投影变换和噪声滤波处理，得到预处理后数据；所述投影变换将地理坐标投影变换为以米为长度单位的平面坐标；
地貌变化差异曲面计算模块包括地貌变化差异曲面确定单元，地貌变化差异曲面确定单元用于将所述预处理后数据D’，以3*3大小的栅格分析窗口为单元，使用空间分析焦点统计功能，统计各单元的标准差值，得到整个地理区域中不同地貌单元标准差值，生成表征地貌高程变化差异程度的栅格化曲面数据Std；
地貌隆起与水平切割深度曲面计算模块包括分水岭提取单元、内接椭圆确定单元、高程算术平均值确定单元、地貌趋势面确定单元、地貌相对高度曲面确定单元和地貌隆起与水平切割深度曲面确定单元；
分水岭提取单元用于对所述预处理后数据D’进行分水岭提取，得到各分水岭地貌单元多边形F_vi，其中i为各分水岭地貌单元的标识号，i＝1……m，m为所述各分水岭地貌单元的总个数；
内接椭圆确定单元用于确定各分水岭地貌单元多边形F_vi的最大内接椭圆Yi和各所述最大内接椭圆的圆心坐标ki、li；
高程算术平均值确定单元用于根据各所述圆心坐标ki、li，确定各所述多边形内的地貌模型栅格点的高程算术平均值ui(ki,li)；
地貌趋势面确定单元用于根据各所述圆心坐标ki、li和所述算术平均值ui(ki,li)采用克里金插值法，得到地貌趋势面W；
地貌相对高度曲面确定单元用于根据所述地貌模型D’和所述地貌趋势面W，将所述地貌模型D’减去所述地貌趋势面W，得到地貌相对高度度量Terr_o；
地貌隆起与水平切割深度曲面确定单元用于将所述地貌相对高度度量Terr_o使用离差标准化法进行归一化处理，得到度量地貌隆起与水平切割深度的栅格化曲面数据Terr；
地形切割密度曲面计算模块包括沟谷线提取单元和地形切割密度曲面确定单元；
沟谷线提取单元用于调入所述预处理后的数据D’，将每个单元海拔值与周边8个单元比较，确定当前单元水流的最大下降方向；累加下降到该点的水文汇流量，储存为F_accum；人机交互确定阈值a，F_accum与a进行大小比较，大于a的单元判断为沟谷地貌，单元值记为1，否则记为Nodata，储存为E，对E矢量化得到整个地理区域上沟谷线数据E_v；根据各所述分水岭地貌单元多边形F_vi与沟谷线数据E_v，进行空间标识，得到表征各分水岭内河流分布的沟谷线数据E_vi,j，其中j为整个地理区域上第i个分水岭地貌单元多边形F_vi内各沟谷线的标识号，j＝1……n，n为所述各分水岭地貌单元内沟谷线的总条数；
地形切割密度曲面确定单元用于根据所述F_vi及E_vi,j数据，进行空间叠加分析和几何特征运算，得到各分水岭地貌单元内沟谷线的总长度Li与所对应的分水岭面积Si；将所述各分水岭地貌单元内沟谷线的总长度Lj除以所述所对应的分水岭面积Si，得到沟壑密度值Mi，公式为Mi＝Li/Si；根据所述整个地理区域上各分水岭地貌单元多边形边界数据F_v和所述沟壑密度值Mi，进行矢量转栅格处理，得到整个地理区域中表征地形切割密度的栅格化曲面数据M。
7.根据权利要求6所述的一种度量地貌破碎程度的栅格化曲面的构建系统，其特征在于：地貌破碎程度曲面计算模块包括地貌破碎程度曲面确定单元，地貌破碎程度曲面确定单元用于将所述地貌高程差异程度的栅格化曲面数据Std、所述地貌隆起与水平切割深度的栅格化曲面数据Terr和所述地形切割密度的栅格化曲面数据M，使用公式A＝Std×Terr×M，进行栅格乘积运算，得到度量地貌破碎程度的栅格化地形曲面A。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110221023.5A CN112950779B (zh) | 2021-02-26 | 2021-02-26 | 一种度量地貌破碎程度的栅格化曲面的构建方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110221023.5A CN112950779B (zh) | 2021-02-26 | 2021-02-26 | 一种度量地貌破碎程度的栅格化曲面的构建方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112950779A CN112950779A (zh) | 2021-06-11 |
| CN112950779B true CN112950779B (zh) | 2021-10-29 |
Family
ID=76246682
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110221023.5A Active CN112950779B (zh) | 2021-02-26 | 2021-02-26 | 一种度量地貌破碎程度的栅格化曲面的构建方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112950779B (zh) |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2003196670A (ja) * | 2001-12-25 | 2003-07-11 | Sony Corp | 侵食地形モデル生成装置、侵食地形モデル生成方法、並びにコンピュータ・プログラム |
| CN108038296A (zh) * | 2017-12-07 | 2018-05-15 | 鲁东大学 | 一种快速估算沟谷沟底线总长度分形计盒维数的方法 |
| CN111854692A (zh) * | 2019-04-26 | 2020-10-30 | 李涛 | 一种在道路测设中的无人机影像匹配点云的测量方法 |
| CN111127646B (zh) * | 2019-12-26 | 2023-03-14 | 西南林业大学 | 一种度量地貌高差的栅格化高程曲面的构建方法及系统 |
-
2021
- 2021-02-26 CN CN202110221023.5A patent/CN112950779B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN112950779A (zh) | 2021-06-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112950777B (zh) | 一种度量地形复杂程度的栅格化曲面的构建方法及系统 | |
| Vivoni et al. | Generation of triangulated irregular networks based on hydrological similarity | |
| CN103363962B (zh) | 一种基于多光谱影像的湖泊水储量遥感估算方法 | |
| CN107180450B (zh) | 一种基于dem的河谷横断面形态的算法 | |
| Gao | Resolution and accuracy of terrain representation by grid DEMs at a micro-scale | |
| Lohani et al. | Application of airborne scanning laser altimetry to the study of tidal channel geomorphology | |
| CN111127646B (zh) | 一种度量地貌高差的栅格化高程曲面的构建方法及系统 | |
| Zhou | Digital elevation model and digital surface model | |
| Werbrouck et al. | Digital elevation model generation for historical landscape analysis based on LiDAR data, a case study in Flanders (Belgium) | |
| Ardiansyah et al. | DEM generation method from contour lines based on the steepest slope segment chain and a monotone interpolation function | |
| CN108648271A (zh) | 一种基于gis数据生成复杂地形网格模型的插值方法 | |
| Szypuła | Digital elevation models in geomorphology | |
| CN104574512A (zh) | 一种顾及地形语义信息的多尺度dem构建方法 | |
| CN117409168B (zh) | 实时动态渲染的洪水预报及洪灾模拟方法和系统 | |
| Šiljeg et al. | The effect of user-defined parameters on DTM accuracy—development of a hybrid model | |
| CN115953556A (zh) | 暴雨内涝道路风险ar预警方法及装置 | |
| Skentos et al. | Landform analysis using terrain attributes. A GIS application on the Island of Ikaria (Aegean Sea, Greece) | |
| Pandey et al. | Urban built-up area assessment of Ranchi township using Cartosat-I stereopairs satellite images | |
| CN112950779B (zh) | 一种度量地貌破碎程度的栅格化曲面的构建方法及系统 | |
| Turel et al. | Delineation of slope profiles from digital elevation models for landslide hazard analysis | |
| Nkeki et al. | Mapping and Geovisualizing Topographical Data Using Geographic Information System (GIS) | |
| Manjare et al. | Digital terrain analysis and geomorphological mapping using remote sensing and GIS: A case study from Central India | |
| CN112950778B (zh) | 一种度量地貌褶皱的栅格化曲面的构建方法及系统 | |
| Gülgen et al. | Automatic extraction of terrain skeleton lines from digital elevation models | |
| Asal | Creation and Analysis of Earth’s Surface Roughness Maps from Airborne LiDAR Measurements in Downtown Urban Landscape |
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 | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20221215 Address after: 518000 211C, Floor 2, West Block, Xincheng Building, No. 1025, Shennan Middle Road, Badeng Community, Nanyuan Street, Futian District, Shenzhen, Guangdong Patentee after: Shenzhen Ring Thumb Electronic Technology Co.,Ltd. Address before: No. 300, bailongsi, Panlong District, Kunming City, Yunnan Province Patentee before: SOUTHWEST FORESTRY University |
|
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20230403 Address after: Zone 13-1-101, Zone A, Longxi No.1, Panlong District, Kunming City, Yunnan Province, 650000 Patentee after: Yunnan outdoor Map Technology Co.,Ltd. Address before: 518000 211C, Floor 2, West Block, Xincheng Building, No. 1025, Shennan Middle Road, Badeng Community, Nanyuan Street, Futian District, Shenzhen, Guangdong Patentee before: Shenzhen Ring Thumb Electronic Technology Co.,Ltd. |