-
中国基本比例尺地形图多采用高斯-克吕格投影(简称高斯投影)。高斯投影为等角横切椭圆柱投影，没有角度变形，中央子午线上没有任何变形，长度比为1，长度比的等变形线平行于中央子午线[1]。由于高斯投影存在面积变形，当精度要求较高时，需要在地球椭球面上进行图斑面积计算，一些文献进行了相关研究[2-7]。这些文献可分为两类：一类是基于坐标直接计算，如文献[8]利用测地坐标计算了椭球面上凸多边形的面积，文献[9]借助空间直角坐标利用曲面积分方法计算了椭球面上的区域面积 ；另一类是基于两条子午线和两条平行圈围成的椭球面梯形间接计算[3, 5, 7, 10-12]。其中，第二次全国土地调查领导小组办公室明确了该类间接计算方法的参数取值。
本文对基于椭球面梯形间接计算图斑面积的方法进行改进，提出了一种不需要梯形面积公式、直接计算图斑面积的方法，并通过实例给出了改进后计算方法的精度判定，以便于进行地球椭球面上的图斑面积计算。
-
文献[10]中基于椭球面梯形间接计算任意图斑面积的方法，分为椭球面梯形面积计算公式和任意图斑椭球面积计算方法两部分，后者大致又可以分为顶层、中间层和底层3层。
顶层给出了整个图斑多边形的计算思路，其本质属于平面多边形面积的计算。对图 1所示多边形的界址点连续编号，计算出其大地坐标P1(B1, L1)、P2(B2, L2)、P3(B3, L3)、P4(B4, L4)，对于任意给定一经线L0，界址点在L0上的投影点为P′1 (B1, L0)、P′2 (B2, L0)、P′3 (B3, L0)、P′4 (B4, L0)，则多边形P1P2P3P4的面积为平面多边形面积计算中常用的4个梯形图块的面积代数和：
$$ \begin{array}{*{20}{c}} {{S_{{P_1}{P_2}{P_3}{P_4}}} = {S_{{P_2}{P_3}{{P'}_3}{{P'}_2}}} + {S_{{P_3}{P_4}{{P'}_4}{{P'}_3}}} + }\\ {{S_{{P_4}{P_1}{{P'}_1}{{P'}_4}}} + {S_{{P_1}{P_2}{{P'}_2}{{P'}_1}}}} \end{array} $$ (1) 
图 1 椭球面上任意图斑的面积计算方法
Figure 1. Area Calculation Method of Patch on Earth Ellipsoid
中间层给出了单个梯形图块椭球面积的计算方法：把整个大梯形图块P1P2P′2P′1拆分为许多个小梯形P1PiP′iP′1，计算其面积Si，Si累加即得到大梯形图块的面积。
底层给出了小梯形椭球面积的计算方法：按椭球面梯形面积计算公式计算，并特意加注参数取值方法，B1、B2分别取沿界址点编号方向的前一个、后一个界址点的大地纬度，ΔL为沿界址点编号方向的前一个、后一个界址点的大地经度的平均值与L0的差。
-
由两条子午线L1、L2和两条纬线B1、B2围成的地球椭球表面称为椭球面梯形(见图 2)。若设经纬度点(B1，L1)附近一个无限小的椭球面小梯形纬度差为dB(单位为rad)，经度差为dL(单位为rad)，对应纬度差弧长为dx(单位为m)，对应经度差弧长为dy(单位为m)，则小梯形面积dF(单位为m2)为：
$$ {\rm{d}}F = {\rm{d}}x \times {\rm{d}}y = M \times N\cos B{\rm{d}}B{\rm{d}}L $$ (2) 
图 2 椭球面梯形
Figure 2. Ellipsoid Trapezoid
式(2)中，M为子午圈曲率半径，N为卯酉圈曲率半径，其计算公式分别为：
$$ M = \frac{{a\left( {1 - {e^2}} \right)}}{{{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{3/2}}}},N = \frac{a}{{{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{1/2}}}} $$ 式中，a为地球椭球体长半径(单位为m)；e为第一偏心率；B为纬度。把M、N的值代入式(2)得：
$$ {\rm{d}}F = \frac{{{a^2}\left( {1 - {e^2}} \right)\cos B}}{{{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^2}}}{\rm{d}}B{\rm{d}}L $$ (3) 把b2=a2(1-e2)代入式(3)，并进行积分，得L1和L2两条经线及B1与B2两条纬线所围成的椭球面梯形面积T：
$$ T = {b^2}\left( {{L_2} - {L_1}} \right)\int_{{B_1}}^{{B_2}} {\frac{{\cos B}}{{{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^2}}}{\rm{d}}B} $$ (4) 对式(4)的被积函数进行幂级数展开，得式(5)，具体推算过程如下。
已知函数f(x)=(1+x)m的幂级数展开式[13]，其中m为任意常数。令x=-(esinB)2，m=-2，并取至x5项，得：
$$ \begin{array}{*{20}{c}} {{{\left( {1 - {e^2}{{\sin }^2}B} \right)}^{ - 2}} = 1 + 2{e^2}{{\sin }^2}B + 3{e^4}{{\sin }^2}B + }\\ {4{e^6}{{\sin }^6}B + 5{e^8}{{\sin }^8}B + 6{e^{10}}{{\sin }^{10}}B} \end{array} $$ 代入式(4)，并分项积分求解得：
$$ \begin{array}{*{20}{c}} {T = {b^2}\left( {{L_2} - {L_1}} \right)\left[ {\sin B + \frac{2}{3}{e^2}{{\sin }^3}B + \frac{3}{5}{e^4}{{\sin }^5}B} \right.}\\ {\left. { + \frac{4}{7}{e^6}{{\sin }^7}B + \frac{5}{9}{e^8}{{\sin }^9}B + \frac{6}{{11}}{e^{10}}{{\sin }^{11}}B} \right]_{{B_1}}^{{B_2}}} \end{array} $$ (5) 对式(5)中正弦函数的高次项sinnB利用三角函数积化和差公式进行降次，得到：
$$ \begin{array}{*{20}{c}} {T = {b^2}\left( {{L_2} - {L_1}} \right)\left[ {A\sin B - B\sin \left( {3B} \right) + C\sin \left( {5B} \right)} \right.}\\ {\left. { - D\sin \left( {7B} \right)E\sin \left( {9B} \right) - F\sin \left( {11B} \right)} \right]_{{B_1}}^{{B_2}}} \end{array} $$ (6) 式(6)中各参数的取值为：
$$ \begin{array}{l} A = 1 + \left( {3/6} \right){e^2} + \left( {30/80} \right){e^4} + \left( {35/112} \right){e^6} + \\ \;\;\;\;\;\;\left( {630/2304} \right){e^8} + \left( {2772/11264} \right){e^{10}} \end{array} $$ $$ \begin{array}{l} B = \left( {1/6} \right){e^2} + \left( {15/80} \right){e^4} + \left( {21/112} \right){e^6} + \\ \;\;\;\;\;\;\left( {420/2304} \right){e^8} + \left( {1980/11264} \right){e^{10}} \end{array} $$ $$ \begin{array}{l} C = \left( {3/80} \right){e^4} + \left( {7/112} \right){e^6} + \\ \;\;\;\;\;\;\left( {180/2304} \right){e^8} + \left( {990/11264} \right){e^{10}} \end{array} $$ $$ D = \left( {1/112} \right){e^6} + \left( {45/2304} \right){e^8} + \left( {330/11264} \right){e^{10}} $$ $$ E = \left( {5/2304} \right){e^8} + \left( {66/11264} \right){e^{10}} $$ $$ F = \left( {6/11264} \right){e^{10}} $$ 把积分上限B2、下限B1代入式(6)，利用三角函数的和差化积公式进行整理，记ΔL=L2-L1，Bn=(B2-B1)/2，Bm=(B2+B1)/2，即可得到式(7)。
$$ \begin{array}{*{20}{c}} {T = 2{b^2}\Delta L\left[ {A\sin {B_n}\cos {B_m} - B\sin \left( {3{B_n}} \right)\cos \left( {3{B_m}} \right)} \right.}\\ { + C\sin \left( {5{B_n}} \right)\cos \left( {5{B_m}} \right) - D\sin \left( {7{B_n}} \right)\cos \left( {7{B_m}} \right) + }\\ {\left. {E\sin \left( {9{B_n}} \right)\cos \left( {9{B_m}} \right) - F\sin \left( {11{B_n}} \right)\cos \left( {11{B_m}} \right)} \right]} \end{array} $$ (7) 综上，式(5)~(7)理论上等价，均为椭球面梯形面积的近似计算公式，计算精度受幂级数展开时保留项数的影响。若忽略e10项，式(7)即为文献[10]中的椭球面梯形面积计算公式；若继续忽略e8项，式(7)即为文献[7]中的椭球面梯形面积计算公式。从原理上讲，两条子午线与两条平行圈所围成的椭球面梯形，无论纬度差大小，都可以直接利用式(7)计算。换言之，对于一个大的椭球面梯形，先沿纬线方向拆分成小的椭球面梯形，再利用式(7)计算，并不会提高面积计算精度。
-
在文献[10]中，其计算方法的中层提出了沿纬线方向切割小梯形问题，并给出了至少需切割为两个部分的最低要求。但是，其对需要切割几个部分以及怎么保证面积计算精度问题等没有给出说明。对图 1中的小梯形进行放大绘图得图 3。
图 3 沿纬线切割的小梯形
Figure 3. Small Trapezoid Cut Along the Weft
图 3中，沿纬线切割得到了小梯形PiPi+1P′i+1P′i，按文献[10]加注的参数取值说明，其实质是用过Pi和Pi+1中点Pm的经线Lm对小梯形进行割、补，把不规则的小梯形面积替换为L0和Lm之间的椭球面梯形面积求解。由于经纬线直交，割、补的两个小直角三角形以相同的纬度差为高，相同的经度差为底，然而，不同纬度、相同经度差的纬线长短不同，两个小三角形面积不等。
小梯形纬度差越大，上、下底边的经度差越大，两个小直角三角形面积差异越大。显然，两个小直角三角形面积的不等是文献[10]要求把整个大梯形沿纬线切割成小梯形进行累加计算的主要原因。
-
在计算任意图斑椭球面积时，为了提高割、补精度，需要尽可能减小切割小梯形的纬度差。当整个大梯形图块切割成的小梯形纬度差很小的时候，用纬差ΔB取代dB，用窄矩形来近似代替窄曲边梯形(小梯形)，从而求得定积分的近似值，即定积分的矩形法近似计算[13]。矩形法定积分近似计算一般先把积分上、下限等分为n个区间，然后再累加计算这n个小梯形的面积。然而，对于不同纬度差的大梯形图块使用相同的n值不利于面积计算精度的控制。本文在近似计算中，使用相同的积分步长ΔB，以适应图斑面积计算中不同大梯形的纬度差差异。无论ΔB取多小，理论上积分上、下限(B2与B1)之间都由整数个小梯形和一个余数小梯形组成，在编程迭代时，循环处理完整数个小梯形后，一定要处理一个余数小梯形。借助定积分的近似计算方法，抛开精度受幂级数展开保留项数影响的椭球面梯形面积公式(7)，可以直接累加计算椭球面梯形面积，化式(4)为式(8)：
$$ T = {b^2}\left( {{L_2} - {L_1}} \right)\sum\limits_{i = {B_1}}^{{B_2}} {\frac{{\cos {B_i}}}{{{{\left( {1 - {e^2}{{\sin }^2}{B_i}} \right)}^2}}}\Delta B} $$ (8) 为适应任意图斑面积的计算，参照文献[10]思路，改造式(8)为式(9)：
$$ T = {b^2}\sum\limits_{i = {B_1}}^{{B_2}} {\frac{{\cos \left( {{B_i} + \Delta B/2} \right)\Delta {L_{i + \Delta B/2}}\Delta B}}{{{{\left( {1 - {e^2}{{\sin }^2}\left( {{B_i} + \Delta B/2} \right)} \right)}^2}}}} $$ (9) 式中，Bi为第i个小梯形的下底边纬度, 其中，第一个小梯形的下底边纬度值为B1；ΔB为整数个小梯形的高(单位为rad)；Bi+ΔB/2为小梯形上、下底边的纬度均值(即图 3中Pm点的纬度值，单位为rad)；ΔLi+ΔB/2是小梯形上、下底边的经度均值Li+ΔB/2与L0子午线经度之差(即图 3中，Pm和P′m之间的纬线长度，单位为rad)。式(9)即为本文所提出的改进算法，因为其使用了图 3中小梯形割、补后的矩形进行近似计算，本文称之为改进的矩形法。
-
不失一般性，本文选取文献[3]上的4个经纬度点(见表 1)围成的椭球面梯形为计算实例，用式(7)计算其椭球面面积, 作为式(8)、式(9)的计算参照。
表 1 西安80坐标系椭球4个点的坐标经纬度[3]
Table 1. Coordinate Latitude and Longitude of Four Points of Xi'an 80 Ellipsoid[3]
点号 1 2 3 4 纬度 39°15′ 39°15′ 39°16′ 39°16′ 经度 116°23′ 116°24′ 116°24′ 116°23′ 中国常见的大地坐标有北京54坐标系、西安80坐标系、2000国家大地坐标系，各大地坐标采用不同的椭球体。在椭球体的长半轴、扁率、地心引力常数和地球自转角速度[14-15]4个常数中，和椭球面面积相关的常数为长半轴a和扁率f。短半轴b和偏心率e为派生参数，计算公式为：b =a×(1-f)；e2=(a2-b2)/a2。为了与原文献计算结果进行比较，本文使用西安80坐标参数：a=6 378 140, f=1/298.257。
在C#语言中，与双精度浮点数相比，十进制数有效数字位数达28~29位[16]，具有更高的精度。本文在累加计算面积过程中，使用十进制数变量以尽可能减少计算误差。由于C#语言提供的三角函数的自变量和返回值为双精度浮点数，在计算时临时把十进制数转换为双精度浮点数调用系统函数，然后再把计算结果转换为十进制数参加其他计算。考虑到三角函数的双精度返回值仅有15~16位有效数字，本文使用“去尾法”保留14位有效数字。
为和椭球面面积进行对照，本文对表 1中的4点，选取117°E为中央经线计算了其高斯投影面积。首先，使用科达普施软件(CODAPS)计算这4个点的高斯坐标，进而计算其投影面积为2 661 913.093 775 950 m2，与文献[3]的计算结果2 661 913.094 m2一致，投影面积比椭球面面积约大180 m2。
-
借助C#语言，采用文献[10]推荐的双精度浮点数变量，式(7)计算的实例梯形的椭球面积为2 661 732.960 119 59 m2；忽略e10项，其椭球面积为2 661 732.960 117 39 m2；继续忽略e8项，其椭球面积为2 661 732.959 433 67 m2；为避免计算失误，同时利用Matlab软件忽略e8项进行了计算复核，MatlaB计算结果为2 661 732.959 432 3 m2。由于三角函数的返回值为双精度浮点数类型，有效数字位数只有15~16位，保留e10项与不保留e10项的计算结果相比，面积精度提高不明显。对照计算结果，判断实例梯形椭球面面积的“去尾法”高精度可靠值为2 661 732.960 11 m2，与文献[3]的计算结果2 661 732.960 m2相一致。国家第二次土地调查要求“面积计算结果以平方米为单位，保留一位小数，四舍五入”，即使忽略e8项，椭球面梯形面积公式仍然满足要求。
取不同ΔB步长，利用式(8)、式(9)分别近似计算实例的椭球面积，其中式(8)分别使用了十进制数类型变量、双精度浮点数类型变量。计算结果见表 2、表 3。
表 2 矩形法定积分近似计算结果
Table 2. Approximate Results of Definite Integral Calculated by Rectangle Method
ΔB取值/rad 椭球面梯形面积/m2 小梯形面积/m2 10-12 2661 732.960 119 3 0.009 151 10-11 2661 732.960 128 9 0.091 514 10-10 2661 732.960 225 2 0.915 143 10-9 2661 732.961 188 4 9.151 434 10-8 2661 732.970 819 9 91.514 34 10-7 2661 733.067 132 1 915.143 4 10-12 2 661 791.843 538 37* 10-10 2 661 732.739 993 85* 10-7 2 661 733.068 534 58* 注：*使用双精度浮点数计算 表 3 改进矩形法定积分近似计算结果
Table 3. Approximation Results of Definite Integral by Improved Rectangle Method
ΔB取值/rad 椭球面梯形面积/m2 小梯形面积/m2 10-10 2 661 732.960 118 2 0.915 143 10-9 2 661 732.960 118 2 9.151 433 10-8 2 661 732.960 118 2 91.514 33 10-7 2 661 732.960 118 2 915.143 3 10-6 2 661 732.960 118 3 9 151.430 10-5 2 661 732.960 129 4 91 513.97 10-4 2 661 732.961 183 6 915 106.6 从理论上讲，由于实例区域位于北半球，积分近似计算方向又是从下限向上限，对于式(8)而言，近似矩形的面积会略微大于其对应的实际椭球面小梯形面积，随着ΔB步长的减少，小梯形面积数值减少，近似计算精度提高。表 2中的计算结果显示，当纬度步长ΔB=10-7 rad时，第一个累加的小梯形面积约915 m2，椭球面梯形面积近似值为2 661 733.067 132 1 m2；当ΔB=10-10时，第一个累加的小梯形面积不足1 m2，椭球面梯形面积近似为2 661 732.960 225 2 m2；当ΔB=10-12时，椭球面梯形面积近似值为2 661 732.960 119 3 m2。随着ΔB的减小，借助十进制数变量，式(8)近似计算结果精度提高，与理论分析相一致。
表 2后3行使用双精度浮点数变量的计算结果显示：当ΔB=10-7时，使用双精度浮点变量的计算精度与使用十进制数变量的计算精度相当；当ΔB=10-10时，随着累加次数千倍的增加，双精度浮点变量有效数字位数不足(仅15~16位)的影响开始显现出来；当ΔB=10-12时，累加误差约60 m2，占高斯投影面积变形误差的1/3。显然，对于积分近似计算的累加公式来说，使用十进制数变量比双精度浮点数变量更适合。
表 3显示，式(9)的计算效率明显优于式(8)。当ΔB=10-4时，式(9)计算的椭球面梯形面积为2 661 732.961 183 6 m2，计算精度与式(8)ΔB=10-9相当；当ΔB=10-5，式(9)计算的椭球面梯形面积为2 661 732.960 129 4 m2，计算精度与式(8)的ΔB=10-11相当；当ΔB=10-7、10-8、10-9、10-10时，使用“去尾法”保留14位有效数字的椭球面梯形面积恒定为2 661 732.960 118 2 m2。将此值与上述ΔB=10-12时式(8)的计算值进行比较，本文判断2 661 732.960 118 2 m2为实例的椭球面积可靠值。
-
如图 1所示，任意图斑多边形P1P2P3P4的面积为4个梯形图块的面积代数和((式1))，其椭球面面积的计算关键在于求得图 1中的大梯形图块的面积，即一条子午线、任意图斑的一条边和两条平行圈围成的梯形面积。
一条子午线(L0=116°22′00″)，表 1中2、4两点的连线，和两条平行圈围成与图 2中P1P2P′2P′1梯形类似的大梯形图块。本文用式(9)计算了该大梯形图块的椭球面面积，具体计算结果见表 4。
表 4 改进矩形法计算的大梯形图块面积
Table 4. Large Ladder Block Area Calculated by Improved Rectangle Method
ΔB取值/rad 小梯形面积/m2 大梯形图块面积/m2 10-11 0.183 028 3 992 651.323 842 9 10-10 1.830 286 3 992 651.323 842 9 10-9 18.302 85 3 992 651.323 842 9 10-8 183.027 1 3 992 651.323 842 8 10-7 1 830.129 3 992 651.323 836 8 10-6 18 287.13 3 992 651.323 230 3 10-5 181 454.9 3 992 651.262 728 9 10-4 1 672 917 3 992 645.527 068 2 表 4显示，随着ΔB的缩小，大梯形图块面积的计算精度提高。当ΔB=10-5时，累加的第一个小梯形面积为181 454.9 m2，大梯形图块面积为3 992 651.262 728 9 m2；当ΔB=10-7时，累加的第一个小梯形面积为1 830.129 m2，大梯形图块面积为3 992 651.323 836 8 m2；当ΔB=10-9时，累加的第一个小梯形面积为1 830.285 m2，大梯形图块面积为3 992 651.323 842 9 m2。当ΔB=10-10、ΔB=10-11时，大梯形图块面积的“去尾法”近似值恒定为3 992 651.323 842 9 m2，该值即为现有C#语言条件下大梯形图块的椭球面积可靠值。
当ΔB=10-4(20.626 5″)时，大梯形图块面积计算值与最终可靠值相差约6 m2，这6 m2误差主要来自于图 3所示两个割、补的小直角三角形面积差异的累加；当ΔB=10-6(0.206 265″)时，大梯形图块面积的割、补误差累计小于0.001 m2。而当大梯形图块不沿纬线方向切割小梯形，而是作为一个小梯形使用，即当ΔB=60″时，大梯形图块面积计算值为3 992 599.440 177 34 m2，比可靠值约小52 m2；当大梯形按文献[10]的最低要求，拆分成两个小梯形，即ΔB=30″时，大梯形图块面积计算值为3 992 638.352 926 54 m2，比可靠值约小13 m2。
综上，式(9)可以回避椭球面梯形的面积计算公式，把任意图斑的面积计算与椭球面梯形的面积计算统一起来，简化任意图斑的面积计算过程。使用十进制数变量，当ΔB=10-6(0.206 265″)时，大梯形图块面积的割、补误差小于0.001 m2；继续缩小ΔB的值，直到小梯形面积的累加值稳定，可以获得大梯形图块的高精度可靠值。
-
1) 椭球面梯形面积公式为近似计算公式，取至e8项，能很好地平衡计算复杂性和计算精度。取至e10项，在现有计算条件下，计算结果精度没有明显提高。
2) 改进的矩形法计算公式(9)式具有很高的计算效率。当ΔB=10-4时，实例中的椭球面梯形面积精度约0.001 m2；当ΔB=10-7时，获得椭球面梯形的高精度可靠值2 661 732.960 118 2 m2。
3) 改进的矩形法计算公式(9)把任意图斑面积计算的中间层和底层统一起来，使得任意图斑的面积计算和椭球面梯形的面积计算共用同一公式，简化了计算过程。使用十进制数变量，缩小ΔB的值，直到小梯形图块的面积累加值稳定，可以获得任意大梯形图块的高精度可靠值，进而可以高精度地计算任意图斑的椭球面面积。
-
摘要: 由于高斯投影存在面积变形，当面积精度要求较高时，需要在地球椭球面上进行图斑面积计算。通过对图斑椭球面积计算方法进行研究，改进了图斑面积的定积分近似计算方法，并利用改进的矩形法进行了实证分析。研究结果显示，与常规矩形法相比，改进的矩形法具有较高的计算效率，可以替代图斑椭球面积计算方法中的中间层算法和底层算法，简化图斑椭球面积的计算过程。借助C#语言的十进制数类型变量，利用改进算法获得了椭球面梯形面积的高精度可靠值。利用改进算法可以计算任意一大梯形图块的椭球面积，进而方便地计算任意图斑的椭球面积。Abstract: Because of the area distortion of Gauss projection, it is necessary to calculate the map patch area on Earth ellipsoid with high area precision in some applications.In this paper, the area calculation method of ellipsoid patch is studied and improved, and the improved algorithm is analyzed empirically. The results show that the approximating definite integrals by the improved Rectangular Rule has higher computational efficiency than by the conventional rectangular rule. The improved rectangle method can replace the middle layer algorithm and the bottom algorithm in ellipsoid patch area calculation method and simplify the calculating process. With the help of decimal number type variables in C # language, a high accuracy value of the ellipsoid trapezoidal area is obtained using the improved algorithm. The ellipsoid area of large trapezoidal block can be obtained, and the ellipsoid area of arbitrary polygon can be easily solved also by the improved algorithm.
-
Key words:
- Earth ellipsoid /
- patch /
- formula for area /
- ellipsoidal trapezoid
-
表 1 西安80坐标系椭球4个点的坐标经纬度[3]
Table 1. Coordinate Latitude and Longitude of Four Points of Xi'an 80 Ellipsoid[3]
点号 1 2 3 4 纬度 39°15′ 39°15′ 39°16′ 39°16′ 经度 116°23′ 116°24′ 116°24′ 116°23′ 
下载: 导出CSV
表 2 矩形法定积分近似计算结果
Table 2. Approximate Results of Definite Integral Calculated by Rectangle Method
ΔB取值/rad 椭球面梯形面积/m2 小梯形面积/m2 10-12 2661 732.960 119 3 0.009 151 10-11 2661 732.960 128 9 0.091 514 10-10 2661 732.960 225 2 0.915 143 10-9 2661 732.961 188 4 9.151 434 10-8 2661 732.970 819 9 91.514 34 10-7 2661 733.067 132 1 915.143 4 10-12 2 661 791.843 538 37* 10-10 2 661 732.739 993 85* 10-7 2 661 733.068 534 58* 注：*使用双精度浮点数计算
下载: 导出CSV
表 3 改进矩形法定积分近似计算结果
Table 3. Approximation Results of Definite Integral by Improved Rectangle Method
ΔB取值/rad 椭球面梯形面积/m2 小梯形面积/m2 10-10 2 661 732.960 118 2 0.915 143 10-9 2 661 732.960 118 2 9.151 433 10-8 2 661 732.960 118 2 91.514 33 10-7 2 661 732.960 118 2 915.143 3 10-6 2 661 732.960 118 3 9 151.430 10-5 2 661 732.960 129 4 91 513.97 10-4 2 661 732.961 183 6 915 106.6
下载: 导出CSV
表 4 改进矩形法计算的大梯形图块面积
Table 4. Large Ladder Block Area Calculated by Improved Rectangle Method
ΔB取值/rad 小梯形面积/m2 大梯形图块面积/m2 10-11 0.183 028 3 992 651.323 842 9 10-10 1.830 286 3 992 651.323 842 9 10-9 18.302 85 3 992 651.323 842 9 10-8 183.027 1 3 992 651.323 842 8 10-7 1 830.129 3 992 651.323 836 8 10-6 18 287.13 3 992 651.323 230 3 10-5 181 454.9 3 992 651.262 728 9 10-4 1 672 917 3 992 645.527 068 2
下载: 导出CSV
-
[1] 黄杏元, 马劲松.地理信息系统概论[M].北京:高等教育出版社, 2008 Huang Xingyuan, Ma Jinsong, An Introduction to Geographic Information System[M].Beijing:Higher Education Press, 2008 [2] 王解先, 俞振武.高斯投影引起的面积计算误差[J].测绘通报, 2003(4):5-6 http://www.cqvip.com/qk/93318X/200304/7732956.html Wang Jiexian, Yu Zhenwu. Area Calculation Error Caused by Gauss Projection[J].Bulletin of Surveying and Mapping, 2003(4):5-6 http://www.cqvip.com/qk/93318X/200304/7732956.html [3] 党亚民, 成英燕, 吴秀娟, 等.不同坐标系图斑理论面积计算研究[J].测绘科学, 2005, 30(6):23-24 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=chkx200506005 Dang Yamin, Cheng Yingyan, Wu Xiujuan, et al. Investigation on the Computation of Theoretic Area of Patch on Map[J]. Science of Surveying and Mapping, 2005, 30(6):23-24 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=chkx200506005 [4] 温珍灵, 崔文刚, 胡君, 等.高斯投影面积变形随偏离中央子午线距离的探讨[J].贵州师范大学学报(自然科学版), 2015, 33(6):16-19 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gzsfdxxb-zr201506003 Wen Zhenling, Cui Wengang, Hu Jun, et al. Gauss Projection Area of Distortion Discussion with Deviation from the Central Meridian[J]. Journal of Guizhou Normal University (Natural Sciences), 2015, 33(6):16-19 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gzsfdxxb-zr201506003 [5] 王玉珍.高斯投影面积变形随纬度变化的研究[J].黑龙江科技信息, 2013, 23:53+55 http://d.wanfangdata.com.cn/Periodical_hljkjxx201323072.aspx Wang Yuzhen. A Study on Area Deformation of Gauss Projection Along with Variation of Latitude[J]. Heilongjiang Science and Technology Information, 2013, 23:53+55 http://d.wanfangdata.com.cn/Periodical_hljkjxx201323072.aspx [6] 刘洋, 高井祥, 王坚, 等.基于抛物线拟合的椭球面区域面积计算方法[J].大地测量与地球动力学, 2015, 35(2):228-231 http://www.cnki.com.cn/Article/CJFDTotal-DKXB201502014.htm Liu Yang, Gao Jingxiang, Wang Jian, et al. The Method of Area Calculation on the Ellipsoid Surface Based on Parabolic Fitting[J].Journal of Geodesy and Geodynamics, 2015, 35(2):228-231 http://www.cnki.com.cn/Article/CJFDTotal-DKXB201502014.htm [7] 张宏, 温永宁, 刘爱利, 等.地理信息系统算法基础[M].北京:科学出版社, 2006 Zhang Hong, Wen Yongning, Liu Aili, et al. Geographic Information System Algorithm Basis[M]. Beijing:Science Press, 2006 [8] 施一民, 朱紫阳.测地坐标计算椭球面上凸多边形面积的算法[J].同济大学学报(自然科学版), 2006, 34(4):504-507 http://www.cnki.com.cn/Article/CJFDTOTAL-TJDZ200604015.htm Shi Yimin, Zhu Ziyang. Computation of Area for Convex Polygon on Earth Ellipsoid by Means of Geodesic Coordinate[J].Journal of Tongji University (Natural Science), 2006, 34(4):504-507 http://www.cnki.com.cn/Article/CJFDTOTAL-TJDZ200604015.htm [9] 林绿, 马劲松.地球椭球面上区域面积的算法研究[J].测绘通报, 2007(6):8-10 http://www.cnki.com.cn/Article/CJFDTOTAL-CHTB200706005.htm Lin Lü, Ma Jinsong. A Method for Area Calculation on Earth Ellipsoid[J]. Bulletin of Surveying and Mapping, 2007(6):8-10 http://www.cnki.com.cn/Article/CJFDTOTAL-CHTB200706005.htm [10] 国务院第二次全国土地调查领导小组办公室. 图幅理论面积与图斑椭球面积计算公式及要求[R]. 国家国土资源部, 北京, 中国, 2008 The Second National Land Survey Leading Group Office of the State Council of the People's Republic of China. The Theoretical Area and the Calculation Formula of Ellipsoidal Area of Patch[R]. Ministry of Land and Resources, Beijing, China, 2008 [11] 程文, 朱磊, 翟宁.基于ArcObject的椭球面图斑面积计算[J].测绘信息与工程, 2010, 35(4):47-48 http://www.cnki.com.cn/Article/CJFDTotal-DKXB201502014.htm Cheng Wen, Zhu Lei, Zhai Ning. Parcel Area Calculation of Ellipsoid Surface Based on ArcObject[J]. Journal of Geomatics, 2010, 35(4):47-48 http://www.cnki.com.cn/Article/CJFDTotal-DKXB201502014.htm [12] 李铁, 王建营.基于不同方法计算椭球面上图斑面积的比较分析[J].经纬天地, 2015(2):18-21 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=xts201502008 Li Tie, Wang Jianying. Comparison and Analysis of Calculated Area of Patch on Ellipsoid Based on Different Methods[J]. Survey World, 2015(2):18-21 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=xts201502008 [13] ESRI. ArcGIS 10. 2. 2 for Desktop Help[R]. ESRI Redlands, California, USA, 2014 [14] 同济大学数学教研室.高等数学[M].北京:高等教育出版社, 1988 Department of Mathematics of Tongji University. Advanced Mathematics[M]. Beijing:Higher Education Press, 1988 [15] 程鹏飞, 文汉江, 成英燕, 等. 2000国家大地坐标系椭球参数与GRS80和WGS84的比较[J].测绘学报, 2009, 38(3):189-194 http://jz.docin.com/p-503223503.html Cheng Pengfei, Wen Hanjiang, Cheng Yingyan, et al. Parameters of the CGCS2000 Ellipsoid and Comparisons with GRS80 and WGS84[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(3):189-194 http://jz.docin.com/p-503223503.html [16] 魏子卿. 2000中国大地坐标系及其与WGS84的比较[J].大地测量与地球动力学, 2008, 28(5):1-5 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dkxbydz200805001 Wei Ziqing. China Geodetic Coordinate System 2000 and its Comparison With WGS84[J]. Journal of Geodesy and Geodynamics, 2008, 28(5):1-5 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=dkxbydz200805001 [17] Microsoft Corporation. Visual Studio (2010-2015) C# Reference[EB/OL]. https://msdn.microsoft.com/en-us/library/678hzkk9.aspx.2016-4-22,2016 -
点击查看大图
图(3) / 表(4)
计量
- 文章访问数: 1478
- HTML全文浏览量: 51
- PDF下载量: 394
- 被引次数: 0
