图斑的椭球面计算研究Word文档下载推荐.docx
《图斑的椭球面计算研究Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《图斑的椭球面计算研究Word文档下载推荐.docx(11页珍藏版)》请在冰豆网上搜索。
(B2﹣B1)—图幅南北图廓的纬差(单位:
弧度),Bm﹦(B1﹢B2)/2。
二、椭球面上任意梯形面积计算公式
(2)
其中:
A,B,C,D,E为常数,按下式计算:
﹢(30/80)е4﹢(35/112)е6﹢(630/2304)е8
﹢(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
米),b—椭球短半轴(单位:
米);
ΔL—图块经差(单位:
弧度);
(B2﹣B1)—图块纬差(单位:
弧度)
Bm﹦(B1﹢B2)/2。
三、高斯投影反解变换()模型
(若坐标不带带号,则不需减去带号×
1000000;
)
+中央子午线经度值(弧度)(3)
公式说明:
若坐标为没有带号前缀格式,则不需减去带号×
若坐标为有带号前缀格式,则需减去带号×
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万
50万
25万
10万
5万
2.5万
1万
5千
经差
6º
3º
1º
30′
15′
7′30″
3′45″
1′52.5″
纬差
4º
2º
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)面积的代数和。
C(B3,L3)
D(B4,L4)
B(B2,L2)
A(B1,L1)
L
L0
C1
D1
A1
Fi
B
Ei(Bi,Li)
图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
ConstZEROAsDouble=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、算法描述
初始化参数
Doublee;
Doublea;
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
计算梯形面积
DoubleAreaVal;
//梯形面积值
DoublelDiference;
//经差
DoublebDiference;
//纬差
DoublebSum;
//纬度和
DoubleItemValue(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*b