《图幅理论面积与图斑椭球面积计算公式及使用 - 中华》由会员分享，可在线阅读，更多相关《图幅理论面积与图斑椭球面积计算公式及使用 - 中华（13页珍藏版）》请在金锄头文库上搜索。
1、图幅理论面积与图斑椭球面积计算公式及使用图幅理论面积与图斑椭球面积计算公式及使用 - - 中华中华自己收集整理的错误在所难免仅供参考交流如有错误请指正！谢谢图幅理论面积与图斑椭球面积计算公式及要求一、 图幅理论面积计算公式（1）式中：a-椭球长半轴(单位：米)-椭球扁率b-椭球短半轴(单位：米)2(a2b2)/a2A1（3/6）2（30/80）4（35/112）6（630/2304）8B （1/6）2（15/80）4（21/112）6（420/2304）8C （3/80）4 （7/112）6（180/2304）8D （1/112）6 （45/2304）8E （5/2304）8L-图幅东西图廓的
2、经差（单位：弧度）(B2B1)-图幅南北图廓的纬差（单位：弧度）Bm（B1B2）/2二、椭球面上任意梯形面积计算公式（2）其中：A,B,C,D,E 为常数按下式计算:2(a2b2)/a2A1（3/6）2（30/80）4（35/112）6（630/2304）8B （1/6）2（15/80）4（21/112）6（420/2304）8C （3/80）4 （7/112）6（180/2304）8D （1/112）6（45/2304）8E （5/2304）8式中：a-椭球长半轴(单位：米)b-椭球短半轴(单位：米)；L-图块经差(单位：弧度)； (B2B1)-图块纬差(单位：弧度)Bm（B1B2）/2三、
3、高斯投影反解变换（）模型（若坐标不带带号则不需减去带号1000000；）+中央子午线经度值（孤度） （3）式中： 公式说明：若坐标为没有带号前缀格式则不需减去带号1000000；若坐标为有带号前缀格式则需减去带号1000000四、计算用到的常数、椭球参数在计算图幅理论面积与任意图斑椭球面积时有关常数及保留的位数按给定数值计算常数：3.14159265358979206264.806247180 椭球常数：= 6378140 = 1/ 298.257= 6356755.29= 6.69438499958795E-03= 6.73950181947292E-03= 6399596.65198801
4、相关常数: k0 = 1.57048687472752E-07k1 = 5.05250559291393E-03k2 = 2.98473350966158E-05k3 = 2.41627215981336E-07k4 = 2.22241909461273E-09五、计算中的取位及要求 高斯投影反解变换后的 B,L 以秒为单位保留到小数点后 6 位四舍五入 采用计算机计算时所有变量数据类型均要定义为双精度 面积计算结果以平方米为单位保留一位小数四舍五入 各种比例尺标准分幅图经差、纬差见表 1 在用大地坐标生成标准分幅图框时要求在每条边框线的整秒处插入加密点表 1 各种比例尺标准分幅图经差、纬差表
5、比例尺1:100 万1:50 万1:25 万1:10 万1:5 万1:2.5 万1:1 万1:5 千经差6o3o1o303015730345152.5纬差4o2o1o20105230115六、任意图斑椭球面积计算方法任意封闭图斑椭球面积计算的原理：将任意封闭图斑高斯平面坐标利用高斯投影反解变换模型将高斯平面坐标换算为相应椭球的大地坐标再利用椭球面上任意梯形图块面积计算模型计算其椭球面积从而得到任意封闭图斑的椭球面积1、计算方法：任意封闭区域总是可以分割成有限个任意小的梯形图块因此任意封闭区域的面积式中 Si 为分割的任意小的梯形图块面积（i=1,2,.n）用公式(2)计算求封闭区域（多边形如图
6、 1）ABCD 的面积 其具体方法为：（1）对封闭区域（多边形）的界址点连续编号（顺时针或逆时针）ABCD提取各界址点的高斯平面坐标 A(X1,Y1)B(X2,Y2)C(X3,Y3)D(X4,Y4)；（2）利用高斯投影反解变换模型公式（3）将高斯平面坐标换算为相应椭球的大地坐标 A(B1,L1)B(B2,L2)C(B3,L3)D(B4,L4)；（3）任意给定一经线 L0（如 L060）这样多边形 ABCD 的各边 AB、BC、CD、DA 与 L0 就围成了 4 个梯形图块（ABB1A1、BCC1B1、CDD1C1、DAA1D1） ； （4）由于在椭球面上同一经差随着纬度升高梯形图块的面积逐渐减
7、小而同一纬差上经差梯形图块的面积相等所以将梯形图块 ABB1A1 按纬差分割成许多个小梯形图块 AEiFiA1用公式（2）计算出各小梯形图块 AEiFiA1 的面积 Si然后累加 Si 就得到梯形图块 ABB1A1 的面积同理依次计算出梯形图块 BCC1B1、CDD1C1、DAA1D1 的面积（注：用公式（2）计算面积时B1、B2 分别取沿界址点编号方向的前一个、后一个界址点的大地纬度L 为沿界址点编号方向的前一个、后一个界址点的大地经度的平均值与 L0 的差） ；（5）多边形 ABCD 的面积就等于 4 个梯形图块（ABB1A1、BCC1B1、CDD1C1、DAA1D1）面积的代数和图 1
8、椭球面上任意多边形计算面积则任意多边形 ABCD 的面积 P 为：P=ABCD= BCC1B1+ CDD1C1+ DAA1D1- ABB1A12、计算要求 利用图形坐标点将高斯坐标系下的几何图形反算投影到大地坐标系进行投影变换 任意指定一条经线 L0从选定多边形几何形状的起始点开始沿顺时针方向依次计算相邻两点构成的线段以及两点到指定经线的平行线构成的梯形面积将该梯形沿纬度变化方向（Y 轴）进行切割至少需切割为 2 个部分 计算过程中应顺同一方向依坐标点逐个计算相邻两点连线与任意经线构成的梯形面积坐标点不得有遗漏若多边形包含内多边形（洞）则该多边形面积为外多边形面积减去所有内多边形面积之和 计算
9、所有梯形面积的代数和即为该多边形的面积七、算法伪代码描述为了确保编程使用的参数、算法一致保证不同软件计算的椭球面积一致我们用算法伪代码描述的方法对编程进行统一在利用计算机编制椭球面积计算软件时计算参数与计算顺序应严格按照以下代码执行1、参数说明双精度类型：圆周率值：PI = 3.14159265358979中央经线:CenterL RHO = 206264.8062471A:ParamAB:ParamBC:ParamCD:ParamDE:ParamEConst ZERO As Double = 0.00000000000180 椭球常数椭球长半轴：aRadius = 6378140 椭球短半轴
10、：bRadius = 6356755.29椭球扁率：ParaAF = 1/ 298.257椭球第一偏心率：ParaE1 = 6.69438499958795E-03椭球第二偏心率：ParaE2 = 6.73950181947292E-03极点子午圈曲率半径：ParaC = 6399596.65198801k0:Parak0 = 1.57048687472752E-07k1:Parak1 = 5.05250559291393E-03k2:Parak2 = 2.98473350966158E-05k3:Parak3 = 2.41627215981336E-07k4:Parak4 = 2.22241
11、909461273E-092、算法描述初始化参数Double e; Double a;e = ParaE2;ParaC = aRadius / (1 - ParaAF);ParamA = 1 + (3 / 6) * e + (30 / 80) * Power(e, 2) + (35 / 112) * Power(e, 3) + (630 / 2304) * Power(e, 4);ParamB = (1 / 6) * e + (15 / 80) * Power(e, 2) + (21 / 112) * Power(e, 3) + (420 / 2304) * Power(e, 4);Param
12、C = (3 / 80) * Power(e, 2) + (7 / 112) * Power(e, 3) + (180 / 2304) * Power(e, 4);ParamD = (1 / 112) * Power(e, 3) + (45 / 2304) * Power(e, 4);ParamE = (5 / 2304) * Power(e, 4);参数初始化结束中央经线转换为弧度CenterL = TransDegreeToArc(CenterL)选定本初子午线为参考经线StandardLat = 0 For 起始点 To 倒数第二点由高斯坐标反解计算经纬度值ComputeXYGeo (P
13、ntColl.Point(i).y, PntColl.Point(i).x, B, L, CenterL)ComputeXYGeo (PntColl.Point(i + 1).y, PntColl.Point(i + 1).x, B1, L1, CenterL)将经纬度转换为弧度值 B = B / RHOL = L / RHOB1 = B1 / RHOL1 = L1 / RHO计算梯形面积 Double AreaVal;/梯形面积值Double lDiference ;/经差Double bDiference; /纬差Double bSum;/纬度和Double ItemValue(5);/计
14、算变量 bDiference = (B1 - B0)；bSum = (B1 + B0) / 2；lDiference = (L1 + L) / 2;ItemValue(0) = ParamA * Sin(bDiference / 2) * Cos(bSum);ItemValue(1) = ParamB * Sin(3 * bDiference / 2) * Cos(3 * bSum);ItemValue(2) = ParamC * Sin(5 * bDiference / 2) * Cos(5 * bSum);ItemValue(3) = ParamD * Sin(7 * bDiference
15、 / 2) * Cos(7 * bSum);ItemValue(4) = ParamE * Sin(9 * bDiference / 2) * Cos(9 * bSum);AreaVal = 2 * bRadius * lDiference * bRadius * (ItemValue(0) - ItemValue(1) + ItemValue(2) - ItemValue(3) + ItemValue(4);areaSum = areaSum + AreaVal; NextEnd Sub3、高斯坐标反解算法Public Sub ComputeXYGeo(x As Double, y As Double, B As Double, L As Double, center As Double)Dim y1 As DoubleDim bf As Doubley1 = y - 500000Dim e As
