基于地理坐标的地理国情统计单元表面面积精确计算

谢成磊，赵荣，梁勇

(1.山东农业大学 信息科学与工程学院，山东 泰安 271018；2.中国测绘科学研究院，北京 100830)

1 引 言

基本统计是地理国情普查工作的重要内容，其根据地理国情普查采集的几何特征类型和地理实体对象，对统计单元内地形地貌、植被覆盖、荒漠与裸露地表、水域、交通网络、居民地与设施等内容进行统计[1]。地表覆盖作为面状要素是地理国情监测的基础，其表面面积是几何特征类型中面状要素在某时刻的一项基本统计指标，不同时间状态下的表面面积能够反映地表覆盖面积的动态变化。目前有很多方法来统计表面积，李轶平等(2007)运用投影面积和坡度的正割值计算表面积[2]，王长海(2010)用逐步分割多变形的方法获得多边形的表面积[3]，江帆(2005)将1个像素格分成2个三角形来统计DEM的表面积[4]，Jeff S.Jenness(2004)用8个周围毗邻点将1个像素格分成8个空间三角形计算DEM的表面积[5]。相比之下，Jeff S.Jenness的工作在处理栅格数据上更具有借鉴意义。这些算法只是提供了基本的计算原理，并没有考虑到地球曲率的影响，且回避了怎样计算目标区域边界附近面积的问题，因此，其并不太适合计算地表面积。本文的算法在前人的基础上做出了适合计算地理表面积的改进。

2 方 法

DEM影像数据是将一系列高程点按规则组织起来的矩阵格网，可用这些高程点来模拟地表的形态特征，计算表面积基本的原理是将DEM组织成三角网，网中的三角形可以近似的代表该三角形区域的面积，所有相邻的三角形组成的三角网的面积可以表示这个DEM的表面积。本文DEM的坐标系统是地理坐标系统。

2.1 原 理

以4个相毗邻的像素格为计算单元，将4个像素的中心点与这4个像素连接起来，组成4个三角形，通过三角形(Ⅰ、Ⅱ、Ⅲ、Ⅳ)的面积计算得到单元的面积。

图1 4个相毗邻的像素格组成一个计算单元示意图

三角形(Ⅰ、Ⅱ、Ⅲ和Ⅳ)的面积可用海伦-秦九韶公式计算。

(1)

其中，a、b、c代表三角形的边长，S是三角形的面积。三角形的边长就是三角形的顶点之间的欧式距离，DEM高程点在地理坐标下，需要转为空间直角坐标才能计算各边长的长度。

图2 过地球旋转轴一个截面的各项参数示意图

图2中O为椭球中心点，a是椭球的长半轴，b为椭球的短半轴，μ为椭圆离心角，θ为椭圆极角，β为地面点Q的纬度，P为过Q点的法线与椭球的交点。

椭圆的参数方程：

(2)

椭圆中各角度参数之间的关系[6]：

(3)

计算地心O到椭球面一点P的距离：

(4)

计算地心O到地表点Q的长，令R=OQ：

(5)

其中，PQ=H=h+ζ，H为大地高，h为正常高，ζ为高程异常。此时可以认为H为DEM的值。

ϑ为地表点Q和地心连线与赤道面的夹角，δ=∠POQ，

(6)

计算地表点Q的空间直角坐标(x，y，z)：

(7)

其中，L为经度。将式(2)到式(6)带入式(7)便可将DEM数据某点的地理坐标转化为空间直角坐标。

空间直角坐标系中任意两点A(xA，yA，zA)和点B(xB，yB，zB)的欧式距离ΔDAB：

(8)

图1中4个栅格值是已知的，而中心点的高程值是不存在的，此时将其转换成空间直角坐标系会发现式(5)、式(6)中的H是未知的，其值需要通过内插得到，可以用反距离加权法[7-8]对中心点的高程值进行估计：

(9)

2.2 流 程

本文的目的是计算多边形圈出范围的表面积，因此在计算时要考虑多边形与栅格点的关系。计算流程如图3所示：

图3 表面积计算的基本流程图

(1)根据多边形的四至范围获取恰好覆盖多边形的DEM数据，DEM的四至范围作为第一个矩形。

(2)对现有的矩形进行分割，矩形的边长长度单位是DEM栅格的长度，在分割时不能分成同样大小的矩形，只能分成4个大小尽可能相同的矩形。

(3)判断所有的矩形与多边形的位置关系。若矩形完全在多边形内，以上文中提到的以单元为单位计算矩形内DEM数据，累加各个计算单元的面积，计算完成后删除矩形；若矩形完全在多边形外，直接删除矩形；若矩形与多边形边界有接触且矩形边长不为1个栅格长度，返回第2步，直至所有的矩形成为边长都为1个像素长度的矩形。

(4)当剩下的矩形的4条边长为一个单位长度时，提取矩形与多边形重合的部分，这部分重新组成小于像素方格的小多边形。对各个小多边形没有高程的边界点用式(9)赋予高程值，用多边形上的点组成带有约束条件的Delaunay三角网[10]，计算并累加每个三角形的面积。得到各个小多边形的面积，累加到多边形内部面积。

计算步骤中的矩形与多边形的位置关系可简单分为三大类：多边形内，多边形外，多边形边界上，根据矩形顶点在多边形内的个数这三类可延伸为6种类型，如图4所示(实心方块代表在多边形内的网格点，灰色部分为多边形)：

图4 矩形与多边形关系示意图

(1)1个网格点在多边形内：矩形只有1个点在多变形内(图4(a))；

(2)2个网格点在多边形内：矩形只有2个相邻点在多边形内(图4(b))；

(3)2个网格点在多边形内：矩形只有2个对角点在多变形内(图4(c))；

(4)3个网格点在多边形内：矩形只有3个点在多边形内(图4(d))；

(5)4个网格点在多边形内：矩形4个点在多边形内，矩形的某条边有可能被多边形边界穿过偶数次，矩形内部可能有多边形边界点(图4(e))；

(6)0个网格点在多边形内：矩形4个点在多边形外，矩形的某条边有可能被多边形边界穿过偶数次，矩形内部可能有多边形的点(图4(f))。

(1)、(2)、(3)和(4)4种情况下矩形与多边形有交集，表示其在多边形的边界上，(5)和(6)情况下矩形有在多边形边界上的可能性，但(5)在一般情况下表示矩形在多变形内，(6)在一般情况下表示矩形在多变形外，但在某些情况下(5)、(6)矩形与多边形是有交集的，此时其在多边形边界上。

矩形与多边形位置关系的判断最基本的是判断矩形顶点与多变形的关系，其基本原理是以影像点为端点向右做出射线，若射线与多边形的交点个数为奇数，则该点在多边形内，否则该点在多边形外。射线法也存在缺点，在使用时需要考虑以下几种特殊情况：

(1)射线恰好穿过多边形最高点或最低点；

(2)射线穿过多边形某条边后，又恰好穿过某个极值点；

(3)射线端点恰好在多边形上；

(4)射线恰好与某一条边重合。

其中，灰色表示点在多边形内，黑色表示点在多边形外。图5右圆形窗口表示某边界处的局部放大图，白线为多边形的边线，黑点为在多变形内部的高程点。

图5 某县县界为对象验证射线法结果图

结合各矩形顶点与多边形的关系以及矩形内是否有多边形的点，便可将矩形分成图4中的6种类型。

3 实验结果

选取中国四川省及附近的20县市为实验区域，分别用本文算法(算法A)和ArcGIS(算法B)计算这些地方的表面面积，结果如表1所示：

表1 算法A与算法B表面积计算结果比较

表1表明算法A所计算的试验区域表面面积的数值平均值比算法B所得数值平均值大约0.61%。

表2表明同一区域DEM分辨率的增加会使该区域的表面积也随之增加。

表3 当图斑尺寸接近DEM分辨率时两种算法计算结果比较

表3和图6表示ArcGIS在处理图斑边界时不精确，当图斑的尺寸达到DEM分辨率的级别时，ArcGIS并不能统计出图斑面积。

ArcGIS计算各县面积需要平面投影后的DEM数据，然后用shp文件中的多边形进行Mask处理，才能统计得到县行政区划的面积。GIS软件通常只有像素方格中心位于多边形内才会将此像素格进行表面积统计，所以Mask处理后的影像边界是锯齿状的，在处理边界问题显得粗糙。用这种方法统计出的地表面积精度和DEM的空间分辨率是息息相关的，DEM分辨率越低，面积统计精度也越低。这种计算方法在处理比较小的图斑的时候会对结果产生较大的影响。图6和表3表示ArcGIS在处理某些图斑边界时不够精确，当图斑的尺寸达到DEM分辨率的级别时，ArcGIS并不能统计出图斑面积。

图6 ArcGIS处理图斑效果图

图6中方块组成的图形是GIS软件根据多边形截取的DEM，其覆盖的范围与多边形有差别。当多边形非常小时，会被GIS软件忽略掉(图6(b)，虚线为1个DEM影像像素的范围，其内部有1个多边形)。

算法A是直接基于地理坐标数据，省去了平面投影过程，避免了投影变形，且算法在边界处完全按照实际边界进行处理，即考虑多边形边界切割DEM影像像素的实际情况进行计算，可以避免投影误差、边界处面积计算误差，并且本算法在图斑面积非常小的情况下也能够计算该图斑的面积。

4 结束语

本文提出的基于地理坐标的表面积算法直接采用大地坐标进行计算，避免了传统软件计算时因平面投影产生的误差传递到表面积统计结果。在处理计算区域边界时，按照边界的实际位置进行处理，使计算结果更加贴近真实值。然而算法也存在不足之处，如目前没有找到更有效率的方法处理多边形内有虫洞的情况，这是以后需要解决的问题。

参考文献：

[1] 国务院第一次全国地理国情普查领导小组办公室.地理国情普查基本统计[M].北京：测绘出版社，2013.

[2] 江帆，吕晓华，王仲兰.基于复化公式的DEM表面积算法分析[J].测绘学院学报，2005，22(4)：263-265.

[3] 王长海.DEM模型的三维地表面积算法研究[J].红水河，2010，29(3)：153-154.

[4] 苗春葆.点与多边形关系的射线法[J].电脑编程技巧与维护，2008(3)：56-58.

[5] 张锦明，郭丽萍，张小丹.反距离加权插值算法中插值参数对DEM插值误差的影响[J].测绘科学技术学报，2012，29(1)：51-56.

[6] JENNESS J S.Calculating landscape surface area from digital elevation models[J].Wildlife Society Bulletin，2004，32(3)：829-839.

[7] AGUILAR F J，AGÜERA F，AGUILAR M A，et al.Effects of terrain morphology，sampling density，and interpolation methods on grid DEM accuracy[J].Photogrmmetric Engineering and Remote Sensing，2005：805-816.

[8] WEINTRIT A.So，What is actually the distance from the equator to the pole?—Overview of the meridian distance approximations[J].The International Journal on Marine Navigation and Safety of Sea Transportation，2013，2(7)：259-272.

[9] TSENG W K，LEE H S.Navigation on a great ellipse[J].Journal of Marine Science and Technology，2010，18(3)：369-375.

[10] SHEWCHUK J R.Delaunay refinement algorithms for triangular mesh generation[J].Computational Geometry，2002，22：21-74.