1、 L图幅东西图廓的经差(单位:弧度)。(B2-Bi)图幅南北图廓的纬差(单位:弧度),Bm= (BmB2) /2o二、椭球面上任意梯形面积计算公式0 1 3S = 2/r AL Asin (B2 B|) cosm - Bsin-(2 -2 2cos3B| + Csin : (2 一 用)cos5Bm7 9-Dsin (Z?2 - B|) cos7+ Esin (Z?2 - ) cos9Z?m(2)(1)其中:A, B, C, D, E为常数,按下式计算:C 2 = (a2 - b2 )/a2A= 1 +(3/6) e 2 + (30/80)e4+ (35/112) e6 +(630/2304)
2、e8(1/6) e2 + (15/80)e4 + (21/112) e6 +(420/2304)(3/80)e4+ (7/112) e6 +(180/2304)(1/112) e 6 +(45/2304)(5/2304)c8米),b一椭球短半轴(单位:米); L图块经差(单位:弧度);(B2-BJ-图块纬差(单位:弧度)Bm= (Bi + B2) /2o三、高斯投影反解变换(模型矿=y-500000-带号x 1000000 (若坐标不带带号,则不需减去带号X 1000000;)E = K*Bf =E + cos E(K sin E-K2 sin3 E + K3 sin5 E-K4 sin7 E
3、)+中央子午线经度值(孤度) (3)t = tgBf 72 = e,2 cos2 Bf N = C/V C = a2/h V + f Ko,Ki,KK3,K4为与椭球常数有关的量。公式说明:若坐标为没有带号前缀格式,则不需减去带号X 1000000;若坐标为有带号 前缀格式,则需减去带号XIOOOOOOo四、计算用到的常数、椭球参数在计算图幅理论面积与任意图斑椭球面积时,有关常数及保留的位数按给定数值计算。常数:JI = 3. 14159265358979p= 206264.806247180椭球常数:轴a = 6378140 椭球扁率 = 1/ 298. 257椭球短半轴= 6356755.
4、 29椭球第二偏心率= 6. 73950181947292E-03极点子午圈川1率半径c = 6399596.65198801相关常数:ko = 1.57048687472752E-07ki = 5.05250559291393E-03k2 = 2.98473350966158E-05k3 = 2.41627215981336E-07k4 = 2.22241909461273E-09五、计算中的取位及要求%1高斯投影反解变换后的B,L以秒为单位,保留到小数点后6位,四舍五入。%1采用计算机计算时,所右变量数据类型均要定义为双精度。%1面积计算结果以平方米为单位,保留一位小数,四舍五入。%1各种
5、比例尺标准分幅图经差、纬差见表1。%1在用大地坐标生成标准分幅图框时,要求在每条边框线的整秒处插入加密点。表1各种比例尺标准分幅图经差、纬差表比例尺1:100 万50 万25 万10 万1 :5万2. 5 万1万5千经差631。30,3015,1! 303, 4552.5”纬差42120,5,2Z 30”15六、任意图斑椭球面积计算方法任意封闭图斑椭球面积计算的原理:将任意封闭图斑高斯平而坐标利用高斯投影反解变 换模型,将高斯平而坐标换算为相应椭球的大地坐标,再利用椭球面上任意梯形图块面积计 算模型计算其椭球血积,从而得到任意封闭图班的椭球血积。1、计算方法:任意封闭区域总是可以分割成有限个任
6、意小的梯形图块,因此,任意封闭区域的而积P = 式中Si为分割的任意小的梯形图块面积(i=l,2,.n)用公式(2)计算。/=1求封闭区域(多边形如图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)任意给定一经线LO (如L0=60。),这样多边形ABCD的各边AB、
7、BC、CD、DA 与 L()就围成了 4 个梯形图块(ABB1A1、BCC1B1、CDD1C1、DAA1D1);(4)由于在椭球面上同一经差随着纬度升高,梯形图块的面积逐渐减小,而同一纬差 上经差梯形图块的面积相等,所以,将梯形图块ABB1A1按纬差分割成许多个小梯形图块 AEiFiAl,用公式(2)计算出各小梯形图块AEiFiAl的面积Si,然后累加Si就得到梯形图 块ABB1A1的面积,同理,依次计算出梯形图块BCC1B1、CDDICk DAA1D1的面积(注: 用公式(2)计算面积时,Bl、B2分别取沿界址点编号方向的前一个、后一个界址点的大 地纬度,AL为沿界址点编号方向的前一个、后一
8、个界址点的大地经度的平均值与L0白笼);(5)多边f? ABCD的面积就等于4个梯形图块(ABB1A1、BCC1B1、CDD1C1、DAA1D1) 而积的代数和。图1椭球面上任意多边形计算面积则任意多边形ABCD的面积P为:P=ABCD= BCC|B|+ CDD|C)+ DAA】D|. ABB)A)2、计算要求%1利用图形坐标点将高斯坐标系下的儿何图形反算投影到大地坐标系,进行投影变换。%1任意指定一条经线L0,从选定多边形儿何形状的起始点开始,沿顺时针方向依次 计算相邻两点构成的线段,以及两点到指定经线的平行线构成的梯形而积。将该梯形沿纬度 变化方向(Y轴)进行切割,至少需切割为2个部分。%
9、1计算过程中应顺同一方向依坐标点逐个计算相邻两点连线与任意经线构成的梯形 面积,坐标点不得有遗漏。若多边形包含内多边形(洞),则该多边形面积为外多边形面积 减去所有内多边形面积之和。%1计算所有梯形面积的代数和即为该多边形的面积。七、算法伪代码描述为了确保编程使用的参数、算法一致,保证不同软件计算的椭球面积一致,我们用算 法伪代码描述的方法对编程进行统一,在利用计算机编制椭球面积计算软件时,计算参数与 计算顺序应严格按照以下代码执行。1、参数说明双精度类型:圆周率值:PI = 3. 14159265358979中央经线:Cent er LRHO = 206264. 8062471A:Param
10、AB:ParamBC:ParamCD:ParamDE:ParamEConst ZERO As Double = 0. 00000000000180椭球常数椭球长半轴:aRadius = 6378140椭球短半轴:bRadius = 6356755. 29椭球扁率:ParaAF = 1/ 298. 257椭球第一偏心率:ParaEl = 6. 69438499958795E-03 椭球第二偏心率:ParaE2 二 6. 73950181947292E-03 极点子午圈曲率半径:ParaC = 6399596.65198801k0:ParakOkl:Parakl=1. 57048687472752
11、E-07=5. 05250559291393E-03k2:Parak2 = 2. 98473350966158E-05k3:Parak3 = 2. 41627215981336E-07k4:Parak4 = 2. 22241909461273E-092、算法描述初始化参数Double e;Double a;e = ParaE2;ParaC = aRadius / (1 - ParaAF);ParamA = 1 + (3 / 6) * o + (30 / 80) * Power (o, 2) + (35 / 112) * Power (e, 3) + (630 / 2304) * Power(e
12、, 4);ParamB = (1 / 6) * e + (15 / 80) * Power (e, 2) + (21 / 112) * Power (e, 3) + (420 / 2304) * Power (c, 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);参数初始化
13、结束中央经线转换为弧度ContcrL = TransDegrecToArc(CcnterL)选定本初子午线为参考经线StandardLat = 0For起始点To倒数第二点由高斯坐标反解计算经纬度值ComputeXYGeo (PntColl. Point(i). y, PntColl. Point(i). x, B, L, CenterL)ComputeXYGeo (PntColl. Point(i + 1).y, PntColl. Point(i + 1). x, Bl, LI, CenterL)将经纬度转换为弧度值B = B / RII0L = L / RHOBl = Bl / RHOLI
14、 = LI / RHO计算梯形面积Double AreaVal;/梯形面积值Double IDiference ; 经差Doubl e bDiference; 纬差Double bSum;/纬度和Double ItcmValue(5) ;/计算变量bDiference =(Bl - BO);bSum = (Bl + BO) / 2;IDiference = (LI + L) / 2;ParamA * Sin (bDi ference / 2) * Cos (bSum):=ParamI3 * Sin (3 * bDi ference / 2) * Cos (3 * bSum);AreaVal =
15、 2 * bRadius * IDiference * bRadius * (TtemValue(0) - TtemValuc(l) + ItemValue(2) - ItcmValuc(3) + ItemValue(4);areaSum 二 areaSum + AreaVal;NextEnd Sub3、高斯坐标反解算法Public Sub ComputeXYGco(x As Double, y As Double, B As Double, L As Double, center As Double)Dim yl As Double Dim bf As Double yl = y - 500
16、000Dim e As Doublee = ParakO * xDim se As Doublese = Sin(e)bf = e + Cos (e) * (Parakl * se - Parak2 * Power (se, 3) + Parak3 * Power (se, 5) - Parak4 * Power (sc, 7)Dim v As DoubleDim t As DoubleDim N As DoubleDim nl As DoubleDim vt As DoubleDim yn As DoubleDim t2 As DoubleDim g As Doublet = Tan (bf
17、)nl = ParaEl * Power (Cos (bf), 2)v = Sqr (1 + nl)N = ParaC / vyn = yl / Nvt = Power (v, 2) * tt2 = 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) / 720B = T ransArcToDegree(B)Dim cbf As Dou
18、blecbf = 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 * 12) * cbf * Power(yn, 5) / 120 + centerL = T ransArcToDegree(L)弧度转换为度Public Function TransArcToDegree(arc As Double) As DoubleDim degree As DoubleDim min As DoubleDim
19、see As DoubleDim ret As DoubleDim tmp As Doubleret = arc 180 / PTdegree = FormatValue(ret, 100, 100) tmp = (ret - degree) * 60min 二 FormtitValue (tmp, 100, 100) sec =(tmp - min) * 60/秒保留到小数点后6位,四舍五入sec = Format(sec, #. 000000)1 EormatValue(sec, 10000000, 100)TransArcToDegree = degree * 3600 + min * 60 + secEnd FunctionPrivate Function FormatValuc(inputVal As Double, procsion As Long, scalcNum As Long) As DoubleFormti t Value 二(Int (inputVcil * precsion) - Int (inputVal * precsion) Mod seal eNum) / precsion
copyright@ 2008-2023 冰点文库 网站版权所有
经营许可证编号:鄂ICP备19020893号-2