[发明专利]基于大地坐标系数字高程模型的地表面积提取方法及系统

权利要求书

1.一种基于大地坐标系数字高程模型的地表表面面积提取方法，其特征在于：进行以下步骤，

步骤一、设原始DEM采用大地坐标系进行存储，DEM栅格像元行数为M，列数为N，栅格分辨率为C，所述DEM为数字高程模型；创建一个(M-1)×2行、(N-1)×2列的空的栅格文件，栅格文件中各栅格像元的值默认为0，栅格文件的坐标系统和原始DEM相同；设原始DEM的起始点经纬度坐标为(B1,L1)，则栅格文件的起始点经纬度坐标为(B1+C,L1+C)；

步骤二、对原始DEM的栅格像元依次遍历，构建规则三角形并计算新的栅格像元区域的表面面积如下，进行以下子步骤，

步骤2.1，设i表示当前遍历的行号、j表示当前遍历的列号，则当前遍历的栅格像元行列坐标为(i,j)，令i初始值为2，j初始值为2，

步骤2.2，根据当前的(i,j)确定当前栅格像元点，并获得行列号为(i-1,j-1),(i,j-1)和(i-1,j)栅格像元点，构成相邻的4个点；

步骤2.3，根据步骤2.2所得相邻的4个点在平面上得到一个四边形，内插出四边形四条边的中点高程值和对角线交点的高程值，得到5个点；

步骤2.4，由步骤2.2所得相邻的4个点和步骤2.3所得4个点相应四边形中四条边的中点和对角线交点，得到空间上相邻的8个三角形；

步骤2.5，对步骤2.2和步骤2.3所得共9个点的经纬度坐标，按照如下公式转换为空间直角坐标，

X＝(N+H)×Cos(B)×Cos(L)

Y＝(N+H)×Cos(B)×Sin(L)

Z＝[N×(1-e²)+H]×Sin(B)

式中，B,L,H分别表示某点的纬度、经度和高程值，X,Y,Z分别表示该点在三维空间直角坐标系中的坐标值；Cos(B)为纬度的余弦值，Sin(B)是纬度的正弦值，Sin(L)是经度的正弦值，Cos(L)是经度的余弦值；e²为地球椭球体的第一偏心率的平方，参数N的计算公式如下，

N=a1-e2sin2B]]>

式中，a为地球椭球体的长半轴；

步骤2.6，利用海伦公式分别计算8个三角形的面积，包括根据三角形三个顶点的空间直角坐标获取面积；

步骤2.7，由相邻三角形两两组合得到四边形，共得到4个新的四边形；

步骤2.8，根据8个三角形与4个新的四边形的组成关系，分别汇总三角形的面积得到4个新的四边形的面积；

步骤2.9，将4个新的四边形的面积写入到新的栅格文件的对应行列号(2i-3,2j-3)、(2i-2,2j-3)、(2i-2,2j-2)、(2i-3,2j-2)中；取原始DEM的下一个遍历的行列坐标为新的(i,j)，返回步骤2.2，直至遍历完成，得到的栅格文件定义为数字地表面积模型，标记为SAM，进入步骤三；

步骤三、给定任意矢量多边形，将矢量多边形与SAM进行空间叠置，汇总SAM中落在矢量多边形内栅格像元的像元值得到V，将V作为该矢量多边形的地表面积；其中，栅格像元是否落在矢量多边形内按照以下规则进行判定，

若栅格像元的中心点落在矢量多边形内或矢量多边形边上，则该栅格像元算在矢量多边形内，否则不在。

2.一种基于大地坐标系数字高程模型的地表表面面积提取系统，其特征在于：包括以下单元，

栅格文件建立单元，用于设原始DEM采用大地坐标系进行存储，DEM栅格像元行数为M，列数为N，栅格分辨率为C，所述DEM为数字高程模型，创建一个(M-1)×2行、(N-1)×2列的空的栅格文件，栅格文件中各栅格像元的值默认为0，栅格文件的坐标系统和原始DEM相同；设原始DEM的起始点经纬度坐标为(B1,L1)，则栅格文件的起始点经纬度坐标为(B1+C,L1+C)；

栅格像元区域表面面积提取单元，用于对原始DEM的栅格像元依次遍历，构建规则三角形并计算栅格像元区域的表面面积，包括以下模块，

初始化模块，用于设i表示当前遍历的行号、j表示当前遍历的列号，则当前遍历的栅格像元行列坐标为(i,j)，令i初始值为2，j初始值为2，

相邻点提取模块，用于根据当前的(i,j)确定当前栅格像元点，并获得行列号为(i-1,j-1),(i,j-1)和(i-1,j)栅格像元点，构成相邻的4个点；

插值模块，用于根据相邻点提取模块所得相邻的4个点在平面上得到一个四边形，内插出四边形四条边的中点高程值和对角线交点的高程值，得到5个点；

三角形构建模块，用于由相邻点提取模块所得相邻的4个点和插值模块所得4个点相应四边形中四条边的中点和对角线交点，得到空间上相邻的8个三角形；

空间直角坐标转换模块，用于对相邻点提取模块和插值模块所得9个点的经纬度坐标，按照如下公式转换为空间直角坐标，

X＝(N+H)×Cos(B)×Cos(L)

Y＝(N+H)×Cos(B)×Sin(L)

Z＝[N×(1-e²)+H]×Sin(B)

式中，B,L,H分别表示某点的纬度、经度和高程值，X,Y,Z分别表示该点在三维空间直角坐标系中的坐标值；Cos(B)为纬度的余弦值，Sin(B)是纬度的正弦值，Sin(L)是经度的正弦值，Cos(L)是经度的余弦值；e²为地球椭球体的第一偏心率的平方，参数N的计算公式如下，

N=a1-e2sin2B]]>

式中，a为地球椭球体的长半轴；

三角形面积提取模块，用于利用海伦公式分别计算8个三角形的面积，包括根据三角形三个顶点的空间直角坐标获取面积；

四边形构建模块，用于由相邻三角形两两组合得到四边形，共得到4个新的四边形；

四边形面积提取模块，用于根据8个三角形与4个新的四边形的组成关系，分别汇总三角形的面积得到4个新的四边形的面积；

数字地表面积模型生成模块，用于将4个新的四边形的面积写入到新的栅格文件的对应行列号(2i-3,2j-3)、(2i-2,2j-3)、(2i-2,2j-2)、(2i-3,2j-2)中；取原始DEM的下一个遍历的行列坐标为新的(i,j)，命令相邻点提取模块工作，直至遍历完成，得到的栅格文件定义为数字地表面积模型，标记为SAM，命令地表面积确定模块工作；

地表面积确定模块，用于给定任意矢量多边形，将矢量多边形与SAM进行空间叠置，汇总SAM中落在矢量多边形内栅格像元的像元值得到V，将V作为该矢量多边形的地表面积；其中，栅格像元是否落在矢量多边形内按照以下规则进行判定，

若栅格像元的中心点落在矢量多边形内或矢量多边形边上，则该栅格像元算在矢量多边形内，否则不在。