图幅理论面积与图斑椭球面积计算公式及要求

一、图幅理论面积计算公式

???-+---??=m

12m 12m 122cos5(2

5Csin cos3(23Bsin cos (21Asin 603604P B B B B B B B B B L πb ）））???-+--m 12m 12cos9(29Esin cos7(27Dsin

B B B B B B ））（1）

式中：

a —椭球长半轴(单位：米)，α—椭球扁率，

b —椭球短半轴(单位：米)。 е2﹦(a 2﹣b 2)/a 2。

A ﹦1﹢（3/6）е2﹢（30/80）е4

﹢（35/112）е6

﹢（630/2304）е8

。 B ﹦ （1/6）е2﹢（15/80）е4

﹢（21/112）е6

﹢（420/2304）е8

。 C ﹦ （3/80）е4

﹢ （7/112）е6

﹢（180/2304）е8

。 D ﹦ （1/112）е6

﹢ （45/2304）е8

。 E ﹦ （5/2304）е8

。 ΔL —图幅东西图廓的经差（单位：弧度）。

(B 2﹣B 1)—图幅南北图廓的纬差（单位：弧度），Bm ﹦（B 1﹢B 2）/2。

二、椭球面上任意梯形面积计算公式

-+---?=m

12m 12m 122cos5(25Csin cos3(23Bsin cos (21Asin 2S B B B B B B B B B L b ）））?

-+--m 12m 12cos9(29Esin cos7(27Dsin B B B B B B ））（2）

其中：A,B,C,D,E 为常数，按下式计算: е2﹦(a 2﹣b 2)/a 2

A ﹦1﹢（3/6）е2﹢（30/80）е4﹢（35/112）е6﹢（630/2304）е8

B ﹦ （1/6）е2﹢（15/80）е4﹢（21/112）е6﹢（420/2304）е8

C ﹦ （3/80）е4﹢ （7/112）е6﹢（180/2304）е8

D ﹦ （1/112）е6﹢（45/2304）е8

E ﹦ （5/2304）е8

式中：a —椭球长半轴(单位：米)，b —椭球短半轴(单位：米)；

ΔL —图块经差(单位：弧度)； (B 2﹣B 1)—图块纬差(单位：弧度) Bm ﹦（B 1﹢B 2）/2。

三、高斯投影反解变换（,,x y B L →）模型

5000001000000y y '=--?带号（若坐标不带带号，则不需减去带号×；） 0E K x =

)sin sin sin sin (cos 7453321E K E K E K E K E E B f -+-+=

()

()()

2424

222222

2'4590617201'935241'21??

? ??++-??? ??-+++??? ??-=N y t V t t N y t V t t N y t V B B f ηη

()

5222223

22'cos 186242851201'cos 12161

'cos 1

??? ?????

??+++++??? ?????

??++-??? ?????

??=N y B t t t N y B t N y B L f f f ηηη +中央子午线经度值（孤度）（3）式中：f tgB t = f B e 222cos '=η V C N /= b a C /2= 21η+=V

。

为与椭球常数有关的量43210,,,,K K K K K 公式说明：若坐标为没有带号前缀格式，则不需减去带号×；若坐标为有带号前缀格式，则需减去带号×。

四、计算用到的常数、椭球参数

在计算图幅理论面积与任意图斑椭球面积时，有关常数及保留的位数按给定数值计算。常数： π﹦3.979

=ρ 206264.

80椭球常数：

a 椭球长半轴 = α椭球扁率 = 1/ 298.257

b 椭球短半轴 = .29

2e 椭球第一偏心率 = 6.95E-03

e '椭球第二偏心率 = 6.92E-03 c 极点子午圈曲率半径 = .

相关常数: k 0 = 1.52E-07 k 1 = 5.93E-03 k 2 = 2.58E-05 k 3 = 2.436E-07 k 4 = 2.73E-09

五、计算中的取位及要求

① 高斯投影反解变换后的B,L 以秒为单位，保留到小数点后6位，四舍五入。 ② 采用计算机计算时，所有变量数据类型均要定义为双精度。 ③ 面积计算结果以平方米为单位，保留一位小数，四舍五入。 ④ 各种比例尺标准分幅图经差、纬差见表1。

⑤ 在用大地坐标生成标准分幅图框时，要求在每条边框线的整秒处插入加密点。

表1 各种比例尺标准分幅图经差、纬差表

六、任意图斑椭球面积计算方法

任意封闭图斑椭球面积计算的原理：将任意封闭图斑高斯平面坐标利用高斯投影反解变换模型，将高斯平面坐标换算为相应椭球的大地坐标，再利用椭球面上任意梯形图块面积计算模型计算其椭球面积，从而得到任意封闭图斑的椭球面积。

1、计算方法：

任意封闭区域总是可以分割成有限个任意小的梯形图块，因此，任意封闭区域的面积

∑==n

i i s P 1

，式中Si 为分割的任意小的梯形图块面积（i=1,2,…n ）用公式(2)计算。

求封闭区域（多边形如图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（如L0＝60°），这样多边形ABCD的各边AB、BC、CD、DA 与L0就围成了4个梯形图块（ABB1A1、BCC1B1、CDD1C1、DAA1D1）；

（4）由于在椭球面上同一经差随着纬度升高，梯形图块的面积逐渐减小，而同一纬差上经差梯形图块的面积相等，所以，将梯形图块ABB1A1按纬差分割成许多个小梯形图块AEiFiA1，用公式（2）计算出各小梯形图块AEiFiA1的面积Si，然后累加Si就得到梯形图块ABB1A1的面积，同理，依次计算出梯形图块BCC1B1、CDD1C1、DAA1D1的面积（注：用公式（2）计算面积时，B1、B2分别取沿界址点编号方向的前一个、后一个界址点的大地纬度，ΔL为沿界址点编号方向的前一个、后一个界址点的大地经度的平均值与L0的差）；

（5）多边形ABCD的面积就等于4个梯形图块（ABB1A1、BCC1B1、CDD1C1、DAA1D1）面积的代数和。

图1 椭球面上任意多边形计算面积

则任意多边形ABCD的面积P为：

P=ABCD= BCC1B1+ CDD1C1+ DAA1D1- ABB1A1

2、计算要求

①利用图形坐标点将高斯坐标系下的几何图形反算投影到大地坐标系，进行投影变换。

②任意指定一条经线L0，从选定多边形几何形状的起始点开始，沿顺时针方向依次计算相邻两点构成的线段，以及两点到指定经线的平行线构成的梯形面积。将该梯形沿纬度变化方向（Y轴）进行切割，至少需切割为2个部分。

③计算过程中应顺同一方向依坐标点逐个计算相邻两点连线与任意经线构成的梯形面积，坐标点不得有遗漏。若多边形包含内多边形（洞），则该多边形面积为外多边形面积减去所有内多边形面积之和。

④计算所有梯形面积的代数和即为该多边形的面积。

七、算法伪代码描述

为了确保编程使用的参数、算法一致，保证不同软件计算的椭球面积一致，我们用算法伪代码描述的方法对编程进行统一，在利用计算机编制椭球面积计算软件时，计算参数与计算顺序应严格按照以下代码执行。

1、参数说明

双精度类型：

圆周率值：PI = 3.979

中央经线:CenterL

RHO = 206264.

A:ParamA

B:ParamB

C:ParamC

D:ParamD

E:ParamE

Const ZERO As Double = 0.

80椭球常数

椭球长半轴：aRadius =

椭球短半轴：bRadius = .29

椭球扁率：ParaAF = 1/ 298.257

椭球第一偏心率：ParaE1 = 6.95E-03

椭球第二偏心率：ParaE2 = 6.92E-03

极点子午圈曲率半径：ParaC = .

k0:Parak0 = 1.52E-07

k1:Parak1 = 5.93E-03

k2:Parak2 = 2.58E-05

k3:Parak3 = 2.436E-07

k4:Parak4 = 2.73E-09

2、算法描述

初始化参数

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);

ParamC = (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 ((i).y, (i).x, B, L, CenterL)

ComputeXYGeo ((i + 1).y, (i + 1).x, B1, L1, CenterL)

将经纬度转换为弧度值

B = B / RHO

L = L / RHO

B1 = B1 / RHO

L1 = L1 / RHO

计算梯形面积

Double AreaVal;//梯形面积值

Double lDiference ;//经差

Double bDiference; //纬差

Double bSum;//纬度和

Double ItemValue(5);//计算变量

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 / 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;

End Sub

3、高斯坐标反解算法

Public Sub ComputeXYGeo(x As Double, y As Double, B As Double, L As Double, center As Double)

Dim y1 As Double

Dim bf As Double

y1 = y - 500000

Dim e As Double

e = Parak0 * x

Dim se As Double

se = Sin(e)

bf = e + Cos(e) * (Parak1 * se - Parak2 * Power(se, 3) + Parak3 * Power(se, 5) - Parak4 * Power(se, 7))

Dim v As Double

Dim t As Double

Dim N As Double

Dim nl As Double

Dim vt As Double

Dim yn As Double

Dim t2 As Double

Dim g As Double

g = 1

t = Tan(bf)

nl = ParaE1 * Power(Cos(bf), 2)

v = Sqr(1 + nl)

N = ParaC / v

yn = y1 / N

vt = Power(v, 2) * t

t2 = Power(t, 2)

B = bf - vt * Power(yn, 2) / 2 + (5 + 3 * t2 + nl - 9 * nl * t2) * vt * Power(yn, 4) / 24 - (61 + 90 * t2 + 45 * Power(t2, 2)) * vt * Power(yn, 6) / 720

B = TransArcToDegree(B)

Dim cbf As Double

cbf = 1 / Cos(bf)

L = cbf * yn - (1 + 2 * t2 + nl) * cbf * Power(yn, 3) / 6 + (5 + 28 * t2 + 24 * Power(t2, 2) + 6 * nl + 8 * nl * t2) * cbf * Power(yn, 5) / 120 + center

L = TransArcToDegree(L)

End Sub

弧度转换为度

Public Function TransArcToDegree(arc As Double) As Double

Dim degree As Double

Dim min As Double

Dim sec As Double

Dim ret As Double

Dim tmp As Double

ret = arc * 180 / PI

degree = FormatValue(ret, 100, 100)

tmp = (ret - degree) * 60

min = FormatValue(tmp, 100, 100)

sec = (tmp - min) * 60

//秒保留到小数点后6位，四舍五入

sec = Format(sec, "####.000000") 'FormatValue(sec, , 100)

TransArcToDegree = degree * 3600 + min * 60 + sec

End Function

Private Function FormatValue(inputVal As Double, precsion As Long, scaleNum As Long) As Double FormatValue = (Int(inputVal * precsion) - Int(inputVal * precsion) Mod scaleNum) / precsion End Function