關于統一圖幅理論面積與圖斑橢球面積計算要求的通知

（國土調查辦發〔2008〕32號）


各省、自治區、直轄市第二次土地調查領導小組辦公室，國土資源廳（國土環境資源廳、國土資源局、國土資源和房屋管理局、房屋土地資源管理局），解放軍土地管理局、新疆生產建設兵團國土資源局：


面積計算是第二次土地調查的一項重要內容，國務院第二次全國土地調查領導小組辦公室組織有關專家，依據《第二次全國土地調查技術規程》，對圖幅理論面積與圖斑橢球面積計算公式進行了細化，明確了面積計算方法，統一了公式中的有關參數，現將《圖幅理論面積與圖斑橢球面積計算公式及要求》予以印發，請各地嚴格遵照執行。



附：圖幅理論面積與圖斑橢球面積計算公式及要求



二〇〇八年三月二十七日

圖幅理論面積與圖斑橢球面積計算公式及要求

一、 圖幅理論面積計算公式
（1）
式中：
a—橢球長半軸(單位：米)，α—橢球扁率，b—橢球短半軸(單位：米)。
е²﹦(a²﹣b²)/a²。
A﹦1﹢（3/6）е²﹢（30/80）е4﹢（35/112）е6﹢（630/2304）е8。
B﹦ （1/6）е²﹢（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—圖幅東西圖廓的經差（單位：弧度）。
(B2﹣B1)—圖幅南北圖廓的緯差（單位：弧度），Bm﹦（B1﹢B2）/2。

二、橢球面上任意梯形面積計算公式
（2）
其中：A,B,C,D,E 為常數，按下式計算:
е²﹦(a²﹣b²)/a²
A﹦1﹢（3/6）е²﹢（30/80）е4﹢（35/112）е6﹢（630/2304）е8
B﹦ （1/6）е²﹢（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—圖塊經差(單位：弧度)； (B2﹣B1)—圖塊緯差(單位：弧度)
Bm﹦（B1﹢B2）/2。

三、高斯投影反解變換（ ）模型
（若坐標不帶帶號，則不需減去帶號×1000000；）



+中央子午線經度值（孤度） （3）
式中：

公式說明：若坐標為沒有帶號前綴格式，則不需減去帶號×1000000；若坐標為有帶號前綴格式，則需減去帶號×1000000。

四、計算用到的常數、橢球參數
在計算圖幅理論面積與任意圖斑橢球面積時，有關常數及保留的位數按給定數值計算。
常數：
π﹦3.14159265358979
206264.8062471
80橢球常數：
= 6378140 = 1/ 298.257
= 6356755.29
= 6.69438499958795E-03
= 6.73950181947292E-03
= 6399596.65198801
相關常數:
k0 = 1.57048687472752E-07
k1 = 5.05250559291393E-03
k2 = 2.98473350966158E-05
k3 = 2.41627215981336E-07
k4 = 2.22241909461273E-09

五、計算中的取位及要求
① 高斯投影反解變換后的B,L以秒為單位，保留到小數點后6位，四舍五入。
② 采用計算機計算時，所有變量數據類型均要定義為雙精度。
③ 面積計算結果以平方米為單位，保留一位小數，四舍五入。
④ 各種比例尺標準分幅圖經差、緯差見表1。
⑤ 在用大地坐標生成標準分幅圖框時，要求在每條邊框線的整秒處插入加密點。
表1 各種比例尺標準分幅圖經差、緯差表
比例尺 1:100萬 1:50萬 1:25萬 1:10萬 1:5萬 1:2.5萬 1:1萬 1:5千
經差 6º 3º 1º30′ 30′ 15′ 7′30″ 3′45″ 1′52.5″
緯差 4º 2º 1º 20′ 10′ 5′ 2′30″ 1′15″

六、任意圖斑橢球面積計算方法
任意封閉圖斑橢球面積計算的原理：將任意封閉圖斑高斯平面坐標利用高斯投影反解變換模型，將高斯平面坐標換算為相應橢球的大地坐標，再利用橢球面上任意梯形圖塊面積計算模型計算其橢球面積，從而得到任意封閉圖斑的橢球面積。
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.14159265358979
中央經線:CenterL
RHO = 206264.8062471

A:ParamA
B:ParamB
C:ParamC
D:ParamD
E:ParamE

Const ZERO As Double = 0.000000000001

80橢球常數

橢球長半軸：aRadius = 6378140
橢球短半軸：bRadius = 6356755.29
橢球扁率：ParaAF = 1/ 298.257
橢球第一偏心率：ParaE1 = 6.69438499958795E-03
橢球第二偏心率：ParaE2 = 6.73950181947292E-03
極點子午圈曲率半徑：ParaC = 6399596.65198801

k0:Parak0 = 1.57048687472752E-07
k1:Parak1 = 5.05250559291393E-03
k2:Parak2 = 2.98473350966158E-05
k3:Parak3 = 2.41627215981336E-07
k4:Parak4 = 2.22241909461273E-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 (PntColl.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 / 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;
Next

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, 10000000, 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