CN109813255A - 一种地球椭球面上大梯形图块的面积计算方法 - Google Patents

Abstract

本发明一种地球椭球面上大梯形图块的面积计算方法，属于地球椭球面任意图斑面积计算领域，本发明通过过大梯形图块斜边中点的经线对大梯形图块进行切割，给出大梯形图块的面积计算公式：椭球面大梯形图块的面积真值等于椭球面梯形面积(近似估计值)减去空白三角形面积(椭球面梯形的缺失部分)，加上阴影三角形面积(从大梯形上割下来的部分)，即，S大梯形图块ABB₁A₁＝S梯形E_AE_BB₁A₁‑SΔEAE_A+SΔEBE_B；而阴影三角形和空白三角形面积可通过递归调用大梯形面积公式完成。当切割三角形面积估计值足够小时，终止递归，获得大梯形图块的高精度面积值。该发明思路清晰，面积精度控制简单可靠。

Description

一种地球椭球面上大梯形图块的面积计算方法

技术领域

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

背景技术

地球椭球面上任意图斑面积计算文献主要分为两类：一类是基于坐标直接计算，如，施一民等(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_iF_iA₁，计算其面积S_i，S_i累加就得到梯形图块ABB₁A₁的面积。底层给出了小梯形图块的计算方法：转换为两条子午线和两条平行圈为界的椭球面梯形计算。两类面积计算方法相比较，后者更易于理解与应用。

第二类方法中，中间层单个大梯形图块的面积计算方法是个难点。史守正等(2018)研究表明，随着拆分小梯形图块的增加，小梯形图块的高越来越小，大梯形图块的面积精度逐渐增加，而对于非常小的小梯形，可以利用改进的矩形直接计算面积，从而回避了椭球面梯形面积公式(史守正，石忆邵，赵伟.椭球面上图斑面积计算方法的改进[J].武汉大学学报·信息科学版，2018，43(5)：779-785)。虽然，这种处理方法合并了中间层算法和底层算法，简化了大梯形图块的椭球面面积计算步骤，但是，也暴露了大梯形拆解为小梯形时需要近似无限拆分的问题，拆分的小梯形的高越小，大梯形图块的面积精度越高。

技术问题：计算任意椭球面图斑面积需要计算大梯形图块面积，而大梯形图块面积计算精度依赖于小梯形的近似无限拆分，这无疑影响了任意图斑椭球面面积计算的应用。目前，国内主流的GIS软件SuperMap及国际知名GIS软件ArcGIS均没有提供可用的图斑椭球面面积的计算功能。

为此，我发明了一种大梯形图块面积计算的新方法，可以简单、高效地计算大梯形图块面积，进而为地球椭球面上任意图斑的面积计算提供支持。

发明内容

不失一般性，以图1中大梯形图块ABB₁A₁为例，说明大梯形图块面积计算的新方法。取AB边的中点为E((B₁+B₂)/2，(L₁+L₂)/2，)，E点沿经线方向在纬线B₁、B₂上的投影点分别为E_A、E_B(见图2)。显然，如果把整个大梯形图块ABB₁A₁当作一个小梯形进行面积计算，那么，计算结果就是椭球面梯形E_AE_BB₁A₁的面积。换句话说，通过经线E_AE_B对原大梯形图块进行切割，用切割掉的经线右边的阴影三角形ΔEBE_B来填补经线左边的空白三角形ΔEAE_A，可以把椭球面大梯形图块转变成为椭球面梯形。切割的阴影三角形与填补的空白三角形都是直角三角形，二者具有相等经度差的底边和相等纬度差的高。然而，由于同等经差低纬度位置的纬线长度大于高纬度位置的纬线长度，空白三角形的面积大于阴影三角形面积，大梯形图块ABB₁A₁的面积实际比椭球面梯形E_AE_BB₁A₁的面积稍小。椭球面大梯形被拆分为若干小梯形，随着拆分小梯形高的减小，切割三角形与填补三角的高减小的同时，他们的底边长度差异也大幅减小，割补后的椭球面梯形面积更接近于椭球面大梯形的面积，从而提高了大梯形图块的面积计算精度。

新方法思路很简单，不走大梯形图块近似拆分为无限多个小梯形的老路，而是直接给出大梯形图块的面积计算真值表达式。大梯形图块面积等于椭球面梯形面积(近似估计值)减去空白三角形面积，再加上阴影三角形面积，即：

S大梯形图块ABB₁A₁＝S梯形E_AE_BB₁A₁-SΔEAE_A+SΔEBE_B

新方法的关键在于对阴影三角形和空白三角形的再审视。切割的阴影三角形和填补的空白三角形均可以视为椭球面大梯形的特例，它们的面积计算方法与常规椭球面大梯形完全一致。阴影三角形和空白三角形的一条经线边可以视为大梯形中的L₀边，二者的一条纬线边对应大梯形的一条纬线边，二者斜边的另一个顶点则对应于大梯形的另一条纬线边，只是该边长度退化为0，这样，割补中产生的阴影三角形和空白三角形就可以被视为短底边退化为0长度的大梯形，因此，也可以通过对斜边的割补，转化为椭球面梯形，使用大梯形图块面积计算真值表达式进行计算(如：图2中的空白三角形进行了进一步割补，取AE边的中点F，过F的子午线在纬线B1上的投影点为F_A，在纬线(B₁+B₂)/2上的投影点为F_E，空白三角形EAE_A被拆分为椭球面梯形EE_AF_AF_E和阴影三角形FAF_A和空白三角形FEF_E)

新方法的实现在于递归算法。在椭球面大梯形面积计算表达式中，椭球面梯形面积计算有现成的计算公式，而两个底、高均缩小的椭球面大梯形面积(阴影三角形和空白三角形作为大梯形的特例，被归类为椭球面大梯形)则通过递归调用原椭球面大梯形图块面积计算表达式进行计算。椭球面梯形面积公式理论上没有误差，随着递归运算，椭球面大梯形面积逐渐减少，大梯形面积近似估计值与真值的误差逐渐减少。当大梯形图块面积估计值小于给定数值时，大梯形面积计算达到给定精度，计算结束，这正符合递归算法的要求。

新方法的计算流程为：

(1)构造函数calTrap(L₁，L₂，B₁，B₂)，计算椭球面梯形面积。该函数使用文献(国务院第二次全国土地调查领导小组办公室.图幅理论面积与图斑椭球面积计算公式及要求[Z].2008-3-28)中的椭球面上任意梯形面积计算公式。具体如下：

其中A、B、C、D、E为常数，按下式计算

A＝1+(3/6)e²+(30/80)e⁴+(35/112)e⁶+(630/2304)e⁸

B＝(1/6)e²+(15/80)e⁴+(21/112)e⁶+(420/2304)e⁸

C＝(3/80)e⁴+(7/112)e⁶+(180/2304)e⁸

D＝(1/112)e⁶+(45/2304)e⁸

E＝(5/2304)e⁸

式中，e为第一偏心率，e²＝(a²-b²)/a²，a为椭球长半轴(单位：m)，b为椭球短半轴(单位：m)，ΔL为梯形图块经差(单位：弧度)，(B₂-B₁)为梯形图块纬差(单位：弧度)，B_m＝(B₁+B₂)/2。

椭球面梯形涉及到两条经线和两条纬线，椭球面梯形面积的正负与该梯形的经度差、纬度差的正负有关。当(L₁-L₂)*(B₁-B₂)＞0时，椭球面梯形获得正的面积值；当(L₁-L₂)*(B₁-B₂)＜0时，椭球面梯形获得负的面积值；当(L₁-L₂)*(B₁-B₂)＝0时，椭球面面积获得0值。

(2)构造递归函数calPatch(L₀，L₁，B₁，L₂，B₂，err)，计算大梯形面积。

与椭球面梯形面积计算中的四个参数不同，椭球面大梯形图块面积计算需要六个参数，以图2为例，L₀为图2中的L₀经线，(B₁，L₁)为A点坐标，(B₂，L₂)为B点坐标，err为精度控制变量，以平方米为单位，比如err＝1m²。

calPatch(L₀，L₁，B₁，L₂，B₂，err)

{

//本函数用来计算椭球面大梯形图块的面积

sum＝calTrap(L₀，(L₁+L₂)/2，B₁，B₂)；//计算椭球面大梯形面积估计值(椭球面梯形面积)

if(Math.Abs(sum)＜err)//当估计值(包含割补的小三角形)的绝对值小于指定值err时，不再递归

return sum；

else//反之，递归计算

return sum+calPatch((L₁+L₂)/2，L₁，B₁，(L₁+L₂)/2，B_m，err)+calPatch((L₁+L₂)/2，(L₁+L₂)/2，B_m，L₂，B₂，err)；//大梯形面积真值＝大梯形面积估计值+阴影三角形面积-空白三角形面积。其中，calPatch((L₁+L₂)/2，L₁，B₁，(L₁+L₂)/2，B_m，err)为空白三角形面积，其值为负，calPatch((L₁+L₂)/2，(L₁+L₂)/2，B_m，L₂，B₂，err)为阴影三角形面积，其值为正。

}

新方法的优势：

相对于把大梯形图块拆解为无数个小梯形图块的旧方法，新方法给出的椭球面大梯形面积公式简单明了；通过精度控制变量err取值和递归算法，新方法可以高效率地适应任意大梯形图块面积的任意精度计算。比如，在特殊情况下，大梯形图块本身为椭球面梯形，旧方法要么需要提前判断，要么仍然继续拆分才能获取高精度面积，而新方法不需要人工干预，只经过一次递归，割、补三角形面积值均为0，就可以获得正确的面积值；对于一般椭球面大梯形，通过更改err变量的值，新方法的面积精度控制更简单，使用该方法，更容易高精度的实现椭球面上任意图斑面积的计算。

附图说明

图1：椭球面上任意多边形计算面积(国务院第二次全国土地调查领导小组办公室.2008)

图2：椭球面上大梯形图块计算面积

图3：椭球面上大梯形图块计算面积实例

具体实施方式

以下将结合具体实施方案来说明本发明。不失一般性，使用西安80椭球参数，选取史守正等(2018)计算过的椭球面大梯形作为实例进行验证。该大梯形两点坐标数值见表1，这两点连线、过两点的两条平行圈和一条子午线(L₀＝116°22′00″)围成了与图2中四边形ABB₁A₁类似的大梯形图块P₂P₄P₅P₆(见图3)。程序代码使用C#语言编制，变量使用十进制变量以尽可能减少计算误差，以史守正等(2018)计算的该大梯形图块面积值3992651.3238429m²作为对照，考察新方法的计算正确性及精度控制变量err对大梯形图块面积精度的控制水平。

计算实例大梯形的面积只需要简单调用递归函数calPatch(L₀，L₁，B₁，L₂，B₂，err)，其中L₀＝116°22′00″，L₁、B₁为点2的经度、纬度，L₂，B₂为点4的经度、纬度，err取值用来控制计算精度，比如err等于1m²。

表1西安80椭球2个点的坐标经纬度

对应不同的面积精度控制变量err的取值，实例大梯形的面积计算值见表2，随着控制变量err数值的减少，面积计算精度逐渐提高，当err取100m²时，实例大梯形的面积可靠值为3992651.32m²时；当err取10m²时，实例大梯形的面积可靠值为3992651.323m²时；当err取1m²时，实例大梯形的面积可靠值为3992651.3238m²时；当err取0.001时，达到实例大梯形的面积可靠值为3992651.3238429m²。该值与史守正等(2018)的计算结果相同，这证明了该方法的正确性。面积精度控制变量err在从100到0.001变化过程中，实例大梯形的面积精度依次提高，err缩小10倍，面积精度提高一位，这种变化规律说明，err变量非常适合椭球面大梯形的面积精度控制。

表2不同精度控制条件下实例大梯形的面积单位：m²

综上，与近似无限切割小梯形的已有算法相比较，新方法思路更简单明了，新方法对大梯形图块的面积精度控制更有效、更方便使用。

Claims (4)

1.一种地球椭球面上大梯形图块的面积计算方法，其特征在于直接给出大梯形图块的面积计算真值表达式：大梯形图块面积等于椭球面梯形面积(近似估计值)减去空白三角形面积，再加上阴影三角形面积。

2.根据权利要求1所述的计算表达式，通过对阴影三角形和空白三角形的再审视，发现椭球面三角形的面积计算可以通过递归调用大梯形图块面积表达式完成。

3.根据权利要求2所述的递归运算，其特征在于通过err变量控制椭球面大梯形面积计算精度。

4.新方法的计算流程为：

(1)构造函数calTrap(L₁，L₂，B₁，B₂)，计算椭球面梯形面积；

(2)构造递归函数calPatch(L₀，L₁，B₁，L₂，B₂，err)，计算大梯形面积，其算法伪代码如下：

calPatch(L₀，L₁，B₁，L₂，B₂，err)

{

//本函数用来计算椭球面大梯形图块的面积

sum＝calTrap(L₀，(L₁+L₂)/2，B₁，B₂)；//计算椭球面大梯形面积估计值(椭球面梯形面积)

if(Math.Abs(sum)＜err)//当估计值(包含割补的小三角形)的绝对值小于指定值err时，不再递归return sum；

else//反之，递归计算

return sum+calPatch((L₁+L₂)/2，L₁，B₁，(L₁+L₂)/2，B_m，err)+calPatch((L₁+L₂)/2，(L₁+L₂)/2，B_m，L₂，B₂，err)；//大梯形面积真值＝大梯形面积估计值+阴影三角形面积-空白三角形面积。其中，calPatch((L₁+L₂)/2，L₁，B₁，(L₁+L₂)/2，B_m，err)为空白三角形面积，其值为负，calPatch((L₁+L₂)/2，(L₁+L₂)/2，B_m，L₂，B₂，err)为阴影三角形面积，其值为正

}。