三调椭球面积计算工具设计与实现

闾海洋，汤 玲，王 刚

(江苏省地质测绘院，江苏 南京 211102)

0 引 言

国务院于2017年10月16日下发《国务院关于开展第三次全国土地调查的通知》(国发〔2017〕48号)，决定自2017年起开展第三次全国土地调查(以下简称“三调”)。为统一三调技术标准，方便各个县区的数据库能够融合到一起，自然资源部于2019年1月28日发布了《第三次全国国土调查技术规程》(TD/T 1055-2019，以下简称《三调技术规程》)。三调以县(区)为单位进行调查，每个县都有一个控制面积，每个县内所有图斑的面积和等于县控制面积，所有县的面积和等于国家的控制面积。通过这种方式进行控制，能够保证全国的面积是一个固定值，不会因为计算方式的不同而产生不同的总面积。三调图斑没有采用高斯投影面上图形的几何面积，而是使用椭球面积作为图斑的面积[1]。计算图斑椭球面积这一功能相对复杂，在主流的GIS建库软件中都没有现成的工具，因此需要使用GIS软件提供的开发接口进行二次开发。

1 面积计算公式解析

《三调技术规程》中对图斑的椭球面积计算方式做出了明确的规定：对任意一个封闭区域(图斑)进行分割，将封闭区域分割成多个梯形图块，先对每个梯形图块进行计算，然后再对所有梯形的面积求和从而得到封闭区域的面积。对于椭球面上的任意一个梯形图块的应采用公式(1)进行计算：

(1)

其中，a表示椭球的长半轴，单位为m;b表示椭球的短半轴，单位为m;ΔL表示梯形图块的经度差值，单位为弧度；BΔ表示梯形图块的纬度差值，单位为弧度；Bm表示梯形图块的南北纬度的平均值，单位为弧度。

《三调技术规程》规定平面坐标使用2000国家大地坐标系，因此参数α、b、e2、A、B、C、D、E都为常数，在计算图斑面积之前可以固化在代码中，具体数值如表1所示。

表1 基本常数表

图斑是土地实际的利用状况在数据库的表现，都是不规则的多边形，因此需要对图斑进行分割处理，将不规则的图斑分割成若干个规则的梯形图块进行计算，分割方式如图1所示。

图1 任意图斑梯形分割示意图

因此图斑ABCDE的总面积公式可以表示为公式(2)：

SABCDE=SABB1A1+SBCC1B1+SCDD1C11+SDEE1D1+SEAA1E1

(2)

式中的每个面积子项都可以通过梯形图块面积计算公式得到。计算梯形图块面积时，使用的坐标是地理坐标，而实际生产过程中使用的坐标是高斯平面坐标，在计算时还需要通过高斯投影反算公式将高斯平面坐标反算为地理坐标。

2 长边问题处理

在计算图斑的椭球面积过程中，为提高面积计算精度，如果图斑节点上任意相邻节点之间的线段长度大于70 m时，需要进行插值计算。内插点仅作为椭球面积计算使用，不需要对线段进行分割。进行内插时，内插点应分布均匀，插值计算流程如图2所示。

图2 插值计算流程

3 面积平差方法

通过面积计算公式计算出每一个图斑的椭球面积，但由于计算过程中会产生一些误差，最终导致面积汇总时总面积会出现不确定性。为解决这一问题，三调采用了控制面积制度。国土调查数据库中的县级及县级以上行政区域界线采用全国陆地行政区域勘界成果确定的界线，乡镇级行政区域界线采用各县(市、区)最新确定的界线。为保证全国各个县(市、区)的国土调查数据库拼到一起后的总面积的确定性，在计算面积时实行了控制面积制度。在开展三调建库时，国家统一下发了各个县(市、区)的控制面积，在一个县(市、区)范围内，所有图斑的总面积应等于下发控制面积。面积计算以平方米为单位，保留2位有效小数。

在三调中，图斑号以村为单位进行编号，控制面积实行三级控制，即县控制各个乡镇的面积，乡镇控制各个村(村级调查区)，村再控制内部的各个图斑。面积平差的总体流程为：① 计算各个乡镇的椭球面积；② 汇总各个乡镇面积，得到汇总面积与县级控制面积的差值，再将差值以乡镇面积大小为权分配到各个乡镇中；③ 计算各个村(村级调查区)的椭球面积；④ 汇总各个村的面积，得到汇总面积与镇控制面积的差值，再将差值以村面积大小为权分配到各个村中；⑤ 计算各个图斑的椭球面积；⑥ 以村为单位进行图斑椭球面积汇总，得到汇总面积与村控制面积的差值，再将差值以图斑面积大小为权分配到各个图斑中。

4 工具设计

椭球面积计算工具分为3个模块来实现，分别为数据管理模块、椭球面积计算模块和面积平差模块。3个模块相互调用，共同实现图斑椭球面积计算的全部功能。

4.1 数据管理模块

三调数据库是一个空间数据库中，行政区(镇)、村级调查区与地类图斑分别保存于XZQ、CJDCQ、DLTB层中，数据结构如表2-表4所示[2]。

表2 XZQ属性结构表(部分)

表3 CJDCQ属性结构表(部分)

表4 DLTB属性结构表(部分)

此模块主要功能有：行政区、村级调查区与地类图斑数据的加载、读取、赋值，要素类属性结构检查，要素的遍历、空间参考获取等操作。使用IWorkspaceFactory.Open方法将三调数据库加载到内存，得到Workspace对象，Workspace对象同时实现了IFeatureWorkspace。使用IFeatureWorkspace. OpenFeatureClass方法代码行政区或地类图斑的名称，得到FeatureClass对象。然后通过IFeatureClass.Search方法获取IFeatureCursor对象后，就可以利用IFeatureCursor.NextFeature方法实现对行政区、地类图斑层中各个要素的遍历[3]。

4.2 椭球面积计算模块

行政区、村级调查区与地类图斑在空间数据库中都是以多边形(Polygon)的形式存在，可以使用统一的方法计算其椭球面积。在GeoDatabase数据模型中，多边形是由一个或多个环(Ring)组成的。环可以分为外环(ExteriorRing)和内环(InteriorRing)，一个外环内可以有0个或多个内环，内环不能独立存在。外环的面积为正数，内环的面积为负数。对于没有岛洞的多边形而言，其只有一个环[4-5]。

椭球面积计算模块主要功能是计算环的椭球面积及多边形内部环的遍历。计算环的椭球面积是此模块中的重点，也是整个工具中最难实现的部分。

对于一个环而言，其主要接口为IRing,但环对象也实现了IPointCollection接口，通过对环中每个节点的遍历，每两个连续点组成一个线段，实现椭球面积的计算，主要实现方法为：① 初始化节点索引N=0 ；② 通过IPointCollection.Point方法，代入索引号得到边 Edeg(N,N+1)；③ 计算Edeg(N,N+1)的长度，如果小于或等于70 m，则直接使用，如果大于70 m，则进行内插，形成一个或多个子边；④ 对于每个子边而言，将首尾结点通过IPoint.Project方法计算点的经纬度坐标；⑤ 通过梯形面积公式(1)计算每个梯形的面积；⑥ 使用N=N+1，回到步骤2再进行计算，直到计算到最后一个点；⑦ 使用公式(2)中提供的方法，汇总所有梯形的面积，得到环的面积。

XZQ层的椭球面积保存于JSMJ字段中，CJDCQ层的椭球面积保存于JSMJ字段中，DLTB层的椭球面积保存于TBMJ字段中。

4.3 面积平差模块

面积平差关系到整个数据库的面积汇总值能否与国家下发的控制面积保持一致。受计算策略、计算精度及小位取位的影响，直接计算得到的面积很难与控制面积保持一致，因此得到图斑的椭球面积后，必须要进行面积平差。进行平差时，分配多出面积的基本策略就以椭球面积值为权进行分配，椭球面积值越大，分配到的误差值就越大。关键算法如下：

(1)汇总XZQ层的JSMJ字段，得到总控制面积与汇总面积的差值，将差值以JSMJ值为权进行分配，平差后的面积保存于DCMJ字段。

(2)一个XZQ内部有一个或多个CJDCQ，以XZQ为单位，分组汇总CJDCQ的JSMJ字段，得到每组的汇总值与相应XZQ的DCMJ值之间的差值，再将差值以CJDCQ的JSMJ值为权进行分配，平差后的面积值保存于CJDCQ的DCMJ字段。

(3)一个CJDCQ内部有一个或多个DLTB，以CJDCQ为单位，分组汇总DLTB的TBMJ，得到每组的汇总值与相应CJDCQ的DCMJ值之间的差值，再将差值以DLTB的TBMJ值为权进行分配，平差后的面积值保存于DLTB的TBMJ字段。

在进行面积平差时，受小数取位的影响，可能一次分配不能达到汇总值与控制值完全相同，还需要进行第二次分配，再次分配还如果还有差值，那么就将这个差值直接分配给面积最大的那个图形。

5 功能实现

椭球面积计算工具的开发基于ArcGIS Add-in框架，IDE选择Visual Studio2013，使用VisualBasic.NET作为编程语言实现以上的全部功能。

所有代码编写、调试完成后，即可对此项目进行编译。编译通过后，会生成一个后缀为esriAddIn的安装包文件。将此文件拷到具有相应ArcGIS版本的电脑上，双击此文件后在弹出的窗口中点击“安装”完成工具部署。打开ArcMap软件，在工具栏中将此工具条显示出来即可使用。图3为椭球面积计算工具的操作界面。

第三次国土调查已经在全国范围内开展，椭球面积计算工具已经应用到了我院承接了多个区县的三调项目中，使用此工具能够快速对县级数据库进行面积计算，计算结果能够通过国家质检软件的检查。

图3 椭球面积计算

6 结 语

本文通过椭球面积计算方法的分析，并以ArcMap AddIn为开发框架，开发了专用的椭球面积计算工具，能够在ArcMap中快速进行椭球面积计算，方便内业人员的建库工作。文中关于ArcMap Add-in开发的案例对开发其他基于ArcMap的应用也有一定的借鉴意义。