利用格网改正法计算椭球面面积

茹仕高，朱紫阳，区永洪，施一民

(1. 广东省国土资源测绘院，广东 广州 510500； 2. 同济大学测绘与地理信息学院，上海 200093)



利用格网改正法计算椭球面面积

茹仕高1，朱紫阳1，区永洪1，施一民2

(1. 广东省国土资源测绘院，广东 广州 510500； 2. 同济大学测绘与地理信息学院，上海 200093)

在对高斯投影面积变形定量分析的基础上，提出了基于高斯投影格网改正的椭球面积计算方法，实现了一次计算、长久受益。该方法直接基于高斯投影面积与已知的格网修正系数，可直接将高斯投影面积转换为椭球面积，计算简便、精度高。分别采用一个大图斑区域和一个小图斑区域进行了计算验证，结果表明，该方法能实现相对精度优于1/100万的计算结果，具有较大的推广价值。

格网改正；椭球面积；投影变形

面积计算是常见的测量工作，高斯投影面积、椭球面积和地表面积是3类常用统计量[1]。高斯投影面积通常的做法是测量边界点的坐标，然后计算封闭区域面积，由各顶点的平面直角坐标计算多边形面积已有简易可行的熟知公式[2]。随着地理信息在各领域的广泛应用，真实地表面与平面差异也需顾及，面积计算由平面转为椭球面，如第二次全国土地调查、第一次全国地理国情普查要求计算图斑椭球面积[3-4]。平面面积与椭球面积相比，投影变形不容忽视，王解先[5]等讨论了高斯投影中央子午线变化和投影面高程变化引起的面积变形及椭球面面积与高斯面上面积的差异；施凤翔[6]分析地籍测量中高斯平面上用解析法计算宗地面积的不足，提出了地籍调查宗地面积也应在椭球面上利用封闭区域界址点的大地坐标精确计算宗地面积的方法；吕长广[7]等分析了土地调查中城镇和农村分别采用投影平面和椭球面计算面积存在的差异，并提出了采用改算系数法进行平面面积改算。然而，按现有方法用大地经纬线计算椭球面多边形面积计算量巨大，为提高计算效率，相关科技工作者进行了大量研究以简化计算量，李玉宝[8]等讨论了高斯投影变形引起的面积测算误差，并分析了变形的量级、性质，提出以高斯平面计算值作为观测值，运用最小二乘原理、拟合出高斯投影变形抛物线函数，对高斯投影面积计算值进行改正；刘洋[9]等在对高斯投影面积变形的性质及量级进行全面分析的基础上，顾及地球曲率，提出抛物线拟合面积修正法，用于计算椭球面区域的面积；施一民[10]等推导了基于大地线与测地平行线正交格网的椭球面三角形面积计算公式。上述研究成果对简化计算量起到了很大的作用，但未就不同抛物线拟合修正适应范围作出分析；测地型椭球面面积主要适用于测地坐标系，对于大地坐标系需要转换，转换计算工作量较大。格网改正法在坐标转换中现已得到广泛的应用[11-13]，那么可否将其引入用于基于高斯投影平面格网改正的椭球面积计算中呢？本文将探索如何采用这种格网改正法来进行简便而又精确的椭球面多边形面积计算。

一、改正原理

1. 高斯投影面积变形

高斯投影为正形投影，根据面积比计算公式可得面积比n与长度比m的关系[14-15]

(1)

式中，Rm为地球平均曲率半径；y为高斯投影横坐标。取至二次项为

(2)

由于高斯投影的正形投影的特性，通过计算高斯投影面上的面积比，则直接可以基于高斯投影面积计算椭球面积。由于高斯投影面为一连续面，计算中不可能也无法实现整个面的面积比计算，也无法通过连续面面积比实现椭球面计算，应用中需进行离散化处理，即格网化处理，在计算误差范围内以格网中心面积比代替格网内任意处面积比。

2. 格网间距确定

从式(2)可以得出，投影平面上面积变形差异最大的方向为沿横轴方向，相近等面积图斑S面积变形差横轴方向上面积都等于S的两图斑的面积变形之差，即

(3)

式中，ym、Δy分别为两相近图斑点高斯投影横坐标均值与差值。为保证计算精度，同名格网内各区域面积变形差应控制在相应的范围υ限内，即υ≤υ限，由于实际应用中取格网中心处面积变形修正值作为格网内区域变形修正值，同名格网内各区域面积变形最大差可放大1倍，即υ≤2υ限，则

(4)

式中，Δy即为满足计算精度要求的最大格网尺寸。在南方地区，若取北回归线(纬度23°26′N)、Rm=6371 km，要求面积变形相对误差均达到百万分之一以下，则测区离开中央子午线的距离与格网间距的关系见表1。

表1　测区离开中央子午线的距离与格网间距的关系

实际应用中，大中比例尺测量与调查采用的是3°带投影，此时测区离开中央子午线的最大距离不超过1.5°，格网间距采用1km即可满足面积变形相对误差达到百万分之一以下。若只需满足面积变形相对误差优于10万分之一，则采用10km的格网间距即可。

3. 格网修正系数计算

(5)

所对应的椭球面上的面积为

(6)

在南方地区，若取北回归线(纬度23°26′N)、2000国家大地坐标系椭球、格网间距为1km，则测区离开中央子午线的距离修正系数的关系见表2。

表2　离开中央子午线的距离与修正系数的关系

由表2可知，格网椭球面上的面积小于高斯投影面积，随着梯形远离中央子午线面积变形值和面积变形相对误差逐渐增大，且与式(2)对应的呈抛物线特性，在高斯投影6°带边缘地区面积变形相对误差达到1/430最大。

同时可知，离投影中央子午线越远，修正系数k值越小。同一格网内，格网中心处修正系数接近真实值；格网中远离中央子午线计算修正系数k值略大于实际修正系数；格网中近中央子午线计算修正系数k值略小于实际修正系数。因此，格网改正法计算格网内图斑面积，离中央子午线近的图斑面积略小于真实图斑面积，离中央子午线远的图斑面积略大于真实图斑面积。

4. 计算步骤

1) 根据计算精度要求，确定格网间距。

2) 根据格网间距，建立覆盖计算区域的格网(若用于软件开发，可计算覆盖整个中国区域)，并计算各格网对应的椭球面上的面积及修正系数。

3) 对图斑与格网进行空间分析，计算图斑对应的格网。

4) 需计算椭球面积多边形的高斯投影面积。

5) 进行计算多边形的面积变形修正，获取各多边形对应的椭球面上的面积。

二、实例分析

1. 算例1

为说明该方法的通用性，选择了如图1所示南方某一县级行政区的遥感调查成果进行图斑椭球面上的面积计算。计算区域总面积约2365km2，离中央子午线最小、最大距离分别为43′、1°25′，共有图斑109 865个，图斑最小上图图斑面积200m2，最大图斑面积235.3km2，平均图斑面积20 712.9m2。按照1km间距建立CGCS2000高斯投影格网，采用10m间距加密格网边界，采用真实椭球面上的面积计算公式分别计算各格网椭球面上的面积及修正系数。

图1　试验区域及图斑边界

为比较计算精度，采用真实椭球面上的面积计算公式计算了各图斑的真实椭球面上的面积，另采用本文介绍的格网改正法计算了各图斑椭球面上的面积，两者比较的统计结果见表3。

表3　图斑格网改正法计算误差统计

计算结果表明，以格网改正法计算的109 865个图斑的椭球面上的面积与理论面积相比较，面积较差绝对值均值为0.017 1m2、相对误差为1/1 210 508，即平均面积变形相对误差在百万分之一以下；最大较差绝对值为33.07m2，其对应相对误差为1/7 114 164，优于700万分之一；最大相对较差为1/175 073，其对应较差绝对值为0.011 7m2，优于10万分之一。

计算结果同时表明，相对精度优于百万分之一的图斑为40 715个，占37.1%，未达到设定的要求百万分之一。究其原因主要是真实椭球面上的面积与改正计算椭球面上的面积差异很小(较差小于1m2的为109 673个图斑，占99.8%；较差大于10m2的图斑数仅为15个)，0.000 1m2的计算误差将产生原本优于百万分之一的相对误差而变为相对误差仅优于10万分之一，而实际绝对较差很小，对面积计算产生的影响完全可以忽略。故整体改正计算效果应是显著的。

2. 算例2

为说明该方法的通用性，选择了广东陆域行政区各县级行政区进行椭球面上的面积计算。计算区域总面积约17.96万km2，选择114°作为高斯投影中央子午线，计算区域离开中央子午线最大距离为4°21′，共有多边形123个，最大多边形面积5634km2，平均多边形面积1 463.13km2。按照10km间距建立CGCS2000高斯投影格网，采用10m间距加密格网边界，采用真实椭球面上的面积计算公式分别计算各格网椭球面上的面积及修正系数。

为比较计算精度，采用真实椭球面上的面积计算公式计算了各图斑的真实椭球面上的面积，并采用本文介绍的格网改正法计算了各图斑椭球面上的面积，得到各图斑的两种面积结果其统计见表4。

表4用格网改正法计算各县级行政区椭球面上的面积的误差统计

类型绝对误差/m2图斑面积/km2相对误差平均较差870.141463.131/1203785最大相对较差1165.37184.651/158444最大绝对较差4193.73857.141/204386

计算结果表明，以格网改正法计算的123个县级行政区的椭球面上的面积与理论面积相比较，面积较差绝对值均值为870.14m2、相对误差1/1 203 785，即平均面积变形相对误差在百万分之一以下；最大较差绝对值为4 193.73m2，其对应相对误差为1/204 386，优于20万分之一；最大相对较差为1/158 444，其对应较差绝对值为1 165.37m2，优于10万分之一。如此大区域计算，都能满足计算精度要求。

计算结果同时表明，相对精度优于百万分之一的图斑为96个，占78.0%；相对精度优于50万分之一的图斑为112个，占91.1%。远优于按照预定的格网间距相对精度应优于10万分之一的要求，究其原因主要有两个方面：①格网间距达到了对相对精度的最低要求；②各图斑跨越格网，格网改正法修正面积与实际椭球面上的面积有正有负，存在计算误差相互抵消的现象。

三、结束语

本文结合地理空间信息应用中大量椭球面上的面积计算的需要，对高斯投影面积变形进行了分析，在此基础上，顾及椭球面上的面积计算的复杂性、高斯投影面积计算的简单性，提出了基于高斯投影格网改正的椭球面上的面积计算方法。该方法直接基于高斯投影面积及已知的格网修正系数，可直接将高斯投影面积转换为椭球面上的面积，计算简便、精度高。

分别采用了一个2000多平方千米的县级行政区遥感调查涉及109 865个图斑、广东省所有县级行政区陆域(含南澳县)123个图斑，分别采用真实椭球面上的面积计算公式与格网改正法计算了相应图斑椭球面上的面积。结果表明，计算结果能实现总体百万分之一的计算精度，能满足应用上的各项需求，该方法具有较大的推广价值。

[1]薛树强，党亚民，秘金钟，等. 顾及非线性地形因子的地表面积计算[J]. 测绘学报，2015，44(3):330-337.

[2]顾孝烈，鲍峰，程效军. 测量学[M]. 4版.上海：同济大学出版社，2011.

[3]中华人民共和国国土资源部.第二次全国土地调查技术规程：TD/T1014—2007[S].北京：[s.n.]，2007.

[4]第一次全国地理国情普查领导小组办公室.地理国情普查基本统计分析技术规定：GDPJ02—2013[S].北京：[s.n.]，2013.

[5]王解先，俞振武. 高斯投影引起的面积误差[J]. 测绘通报，2003(4)：5-6.

[6]施凤翔. 基于椭球面的城镇地籍宗地面积量算方法研究[C]∥地理信息与物联网论坛暨江苏省测绘学会2010年学术年会论文集. 无锡：江苏省测绘学会，2010.

[7]吕长广，张先为，李淑青. 土地调查中城乡土地面积计算的差别及其处理方法[J]. 地矿测绘，2008，24(3)：9-11.

[8]李玉宝，曹智翔. 顾及地球曲率的面积计算问题[J]. 重庆交通学院学报，2007，26(1)：152-154.

[9]刘洋,高井祥,王坚，等. 基于抛物线拟合的椭球面区域面积计算方法[J]. 大地测量与地球动力学, 2015, 35(2): 228-231.

[10]施一民，朱紫阳. 测地坐标计算椭球面上凸多边形面积的算法[J]. 同济大学学报(自然科学版)，2006,34(4):504-507.

[11]郭充，吕志平，李岩，等. 基于格网的坐标转换方法[J]. 信息工程大学学报，2010，11(2)：166-169.

[12]吕志平，魏子卿，李军，等．CGCS2000高精度坐标转换格网模型的建立[J]. 测绘学报，2013，42(6):791-797.

[13]刘文建，施闯. 基于格网改正的CORS实时三维坐标服务系统[J]. 测绘通报，2013(10)：70-72.

[14]施一民. 现代大地控制测量[M]. 2版.北京：测绘出版社，2007.

[15]孔祥元，郭际明，刘宗泉. 大地测量学基础[M]. 2版.武汉：武汉大学出版社，2010.

Study on Ellipsoid Surface Area Calculation Method Using Grid Correction Technology

RU Shigao，ZHU Ziyang，OU Yonghong，SHI Yimin

10.13474/j.cnki.11-2246.2016.0274.

2015-10-08；

2016-03-07

国家科技支撑计划(2012BAB16B01)；国家自然科学基金(40471114)

茹仕高(1965—)，男，高级工程师,主要研究方向为大地测量和工程测量。E-mail:332575991@qq.com

P226

B

0494-0911(2016)08-0124-04

引文格式：茹仕高，朱紫阳，区永洪，等.利用格网改正法计算椭球面面积[J].测绘通报，2016(8)：124-127.