[发明专利]一种地球椭球面上任意图斑的面积计算方法

说明书

本发明一种地球椭球面上任意图斑的面积计算方法，算法思路分顶层、中间层和底层：底层本质上是椭球面梯形的面积计算，本方法借助十进制数变量及函数，打破椭球面梯形面积精确计算公式计算精度局限，获得高达23位有效数字的面积值，保证了任意图斑的面积计算精度；顶层分解任意图斑为各椭球面大梯形，通过各大梯形面积的代数和累加获得任意图斑面积；中间层则通过递归调用底层算法获取大梯形图块的面积。本方法给出了任意图斑面积计算的完整流程，并通过实例进行了验证，使顺时针多边形获得正面积，逆时针多边形获得负面积。本方法计算思路明了，图斑面积精度控制简单可靠，易于软件实现及推广。

技术领域

本发明属于地球椭球面任意图斑面积计算领域，具体地说，为了避免高斯-克吕格投影变形的影响，在面积精度要求较高的应用中，需要计算土地斑块的地球椭球面面积时，使用该发明方法。

背景技术

地球椭球面上任意图斑面积计算文献主要分为两类：一类是基于坐标直接计算，如，施一民等(2006)利用测地坐标计算了椭球面上凸多边形的面积(施一民，朱紫阳.测地坐标计算椭球面上凸多边形面积的算法[J].同济大学学报(自然科学版)，2006，34(04)：504-507)，林绿等(2007)借助空间直角坐标利用曲面积分的方法对椭球面上区域面积的计算方法进行了研究(林绿，马劲松.地球椭球面上区域面积的算法研究[J].测绘通报，2007，(06)：8-10)；另一类是基于椭球面上两条子午线和两条平行圈为界的梯形间接计算(国务院第二次全国土地调查领导小组办公室.图幅理论面积与图斑椭球面积计算公式及要求[Z].2008-3-28)，该计算方法大致可以分为三个层次：顶层、中间层和底层。顶层给出整个图斑多边形的计算思路，即：多边形ABCD的每一条边与任意给定经线L₀围成一个大梯形图块(见图1，以AB边为例，A、B点两点沿纬线方向在经线L₀上的投影点分别为A₁、B₁，则四边形ABB₁A₁即为AB边与L₀围成的大梯形图块)，对围成的所有大梯形图块面积计算其代数和，就得到多边形ABCD的面积。中间层给出了单个大梯形图块的计算方法：把单个大梯形图块ABB₁A₁按纬线切割，拆分成许多个小梯形图块AE₁F₁A₁，计算其面积S₁，S₁累加就得到梯形图块ABB₁A₁的面积。底层给出了小梯形图块的计算方法：转换为两条子午线和两条平行圈为界的椭球面梯形计算，使用近似计算公式计算其面积。两类面积计算方法相比较，后者更易于理解与应用。然而，第二类方法，并没有给出大梯形具体拆分的细节，具体应该拆分成多少个小梯形？或者拆分后的小梯形高等于多少合适？史守正等(2018)对该问题进行了探索性研究，研究表明，随着拆分小梯形图块的高越来越小，可以采用定积分近似计算的思路替代近似计算公式进行大梯形图块的面积计算，且计算精度逐渐增加，并获得样本梯形的高精度可靠值2661732.9601161m²(史守正，石忆邵，赵伟.椭球面上图斑面积计算方法的改进[J].武汉大学学报·信息科学版，2018，43(5)：779-785)。

技术问题：计算任意椭球面图斑面积需要计算大梯形图块面积，而大梯形图块面积计算精度依赖于小梯形的近似无限拆分，以及小梯形的面积计算。椭球面小梯形面积计算实质上是利用椭球面梯形来完成的，椭球面梯形面积公式有精确计算公式和近似计算公式之分，在使用自变量和返回值均为双精度浮点数的系统函数的前提下，采用近似计算公式计算精度较高。然而，双精度浮点数只有15-16位有效数字，对于单个椭球面梯形面积计算而言精度尚可，对于近似无限拆分的小梯形而言，该精度可能影响其累加后的大梯形面积精度。迄今为止，椭球面上任意图斑的面积计算仅仅局限于理论探索阶段，国内主流的GIS软件SuperMap及国际知名GIS软件ArcGIS均没有提供可用的图斑椭球面面积的计算功能。