测量专业基本计算与编程

在赤道时算法不能正常计算。电脑的广泛应用为改进大地主题的解算方法提供了新的条件，在电算中进行迭代计算，简单易行，这里介绍一种嵌套系数法精密解算任意距离的大地主题问题，仅作为教材内容的补充。 1、符号意义

， 参考椭球的长、短半轴； 参考椭球的扁率； 参考椭球的二偏心率；

、 椭球面上两点的大地纬度； 、 椭球面上两点的大地经度； 两点大地经差；

椭球面上两点的大地线长度； 、 两点的大地线方位角； 、 两点的归化纬度；

大地线S 的最高纬度； 两点的球面经差； 两点的球面角距；

大地线自赤道交点至1点间的球面角距； 大地线自赤道交点至2点间的球面角距；

大地线自赤道交点至1、2点间的球面角距之和 ； 球面经差 的改正项； 球面经差 的改正项； 、 计算嵌套系数的算子；

、 计算椭球面距离改正项的嵌套系数；

计算椭球面经差改正项的嵌套系数。 2、嵌套系数法大地主题正算 已知： ，计算 。 1）计算准备

2）嵌套系数计算 ，

3）球面角距 改正项 的迭代计算； （首次取 ） 或 ；

采取迭代算法，直到 的变化小于限差为止，一般取限差为 弧度。 4）计算 ； 或

，在此需要判断象限； ； ；

，在此需要判断象限； 5） 检验计算 。

3、嵌套系数法大地主题反算 已知： ，计算 。 1）计算准备

2）球面经差 的改正数 的计算 （首次取 ） ，并判断象限 ；

采取迭代算法，直到 的变化小于限差为止，一般取限差为 弧度。 3）球面距离 改正项 的迭代计算

4）计算 、 、 ；

，并判断象限； ，并判断象限； 5）检核计算 。

4、嵌套系数法大地主题算例与计算程序 1）正算算例

大地主题正算

================== 计算方法：嵌套系数法( 适合任意距离 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 30.29582043 120.05402184 247.27504300 28230.9350 30.24058365 119.49233864 67.19353763

3 4 47.46526740 35.49363300 44.12136600 44797.2820 48.04096664 36.14450516 224.30535480

5 6 40.02356784 130.10122676 1.49430000 80000.0000 40.45479037 130.12011050 181.50535462

7 8 40.02356784 115.10000000 36.12010300 414306.5380 43.00558793 118.10030013 218.11268010

9 10 68.58000000 33.05000000 339.49563900 7999606.3800 37.44599774 -122.26000117 9.01078071

11 12 35.00002200 90.00001100 100.00003300 15000000.20 -30.29209655 215.59043332 290.32533894

2）反算算例

大地主题反算

================== 计算方法： 嵌套系数法( 适合任意距离 ) 参考椭球：克拉索夫斯基椭球

N01 NO2 B1/B2 L1/L2 A12/A21 S12 1 2 30.29582043 120.05402184 247.27504216 28230.9375 30.24058354 119.49233853 67.19353669

3 4 47.46526470 35.49363300 44.12102618 44796.5558 48.04096384 36.14450004 224.30501114

5 6 40.02356784 130.10122627 1.49432980 80000.0045 40.45479027 130.12011040 181.50538463

7 8 40.02356784 115.10000000 36.12010267 414306.5361 43.00558784 118.10030000 218.11267964

9 10 68.58000000 33.05000000 339.49563867 7999606.3881 37.44599755 -122.26000057 9.01078088

11 12 35.00002200 90.00001100 100.00003279 15000000.3385 -30.29209640 215.59043380 290.32533886

注：上述正反算数据例1取自朱华统：《椭球大地计算》P100算例；例2自陈键、晁定波主编《椭球大地测量学》P98面的数据；例3、例4取自赵文光著《椭球大地测量学》P88和P93例，例5取自陈键、晁定波主编《椭球大地测量学》P111和P113面的数据。例54取自周江华：《测绘通报》，2002，5（6）：108－111。

3）嵌套系数法正反算计算程序

Sub DDZT_Nest(B1 As Double, L1 As Double, A1 As Double, S12 As Double, _ B2 As Double, L2 As Double, A2 As Double, k) eps = e2 / (1# - e2) '计算椭球二偏心率 alpha = 1 - Sqr(1 - e2) b = a / Sqr(1# + eps) Select Case k Case 1 '正算

u1 = Atn(Sqr(1 - e2) * Tan(B1)) sigma1 = Atn(Tan(u1) / Cos(A1))

cosUn = Cos(u1) * Sin(A1): sinUn2 = 1 - cosUn * cosUn

t = eps * sinUn2 / 4 V = alpha * sinUn2 / 4

xk1 = 1 + t * (1 - t * (3 - 5 * t + 11 * t * t) / 4) xk2 = t * (1 - t * (2 - t * (37 - 94 * t) / 8))

xk3 = V * (1 + alpha + alpha * alpha - V * (3 + 7 * alpha - 13 * V)) Dsigma = 0 Do

Dsigma0 = Dsigma

sigma = S12 / xk1 / b + Dsigma0 Msigma = 2 * sigma1 + sigma

xx = Cos(sigma) * Cos(2 * Msigma) + xk2 * (1 + 2 * Cos(2 * sigma)) _ * Cos(3 * Msigma) / 6

Dsigma = xk2 * Sin(sigma) * (Cos(Msigma) + xk2 * xx / 4) Loop While Abs(Dsigma - Dsigma0) * P0 > 0.00001

sinB2 = Sin(u1) * Cos(sigma) + Cos(u1) * Sin(sigma) * Cos(A1) sigma2 = sigma1 + sigma

cosB2 = Sqr(1 - e2) * Sqr(1 - sinUn2 * Sin(sigma2) * Sin(sigma2)) B2 = Atn(sinB2 / cosB2)

'sinU2 = Sin(U1) * Cos(sigma) + Cos(U1) * Cos(A1) * Sin(sigma) 'B2 = Atn(sinU2 / Sqr(1 - e2) / Sqr(1 - sinU2 * sinU2)) sinw = Sin(sigma) * Sin(A1)

cosw = Cos(u1) * Cos(sigma) - Sin(u1) * Sin(sigma) * Cos(A1) tanw = sinw / cosw: w = Abs(Atn(sinw / cosw)) sinA1 = Sin(A1)

If tanw > 0 And sinA1 > 0 Then w = w If tanw < 0 And sinA1 > 0 Then w = pi - w If tanw < 0 And sinA1 < 0 Then w = -w If tanw > 0 And sinA1 < 0 Then w = w - pi

dw0 = sigma + xk3 * Sin(sigma) * (Cos(Msigma) + xk3 * Cos(sigma) * Cos(2 * Msigma))

dw = (1 - xk3) * alpha * cosUn * dw0 L2 = L1 + w - dw

SinA2 = Cos(u1) * Sin(A1)

cosA2 = Cos(u1) * Cos(sigma) * Cos(A1) - Sin(u1) * Sin(sigma) tanA2 = SinA2 / cosA2: A2 = Abs(Atn(SinA2 / cosA2)) If sinA1 < 0 And tanA2 > 0 Then A2 = A2 If sinA1 < 0 And tanA2 < 0 Then A2 = pi - A2 If sinA1 > 0 And tanA2 > 0 Then A2 = pi + A2 If sinA1 > 0 And tanA2 < 0 Then A2 = 2 * pi - A2 H = Cos(u1) * sinA1 - Cos(u2) * Sin(A2) '检核计算 Case 2 '反算

u1 = Atn(Sqr(1 - e2) * Tan(B1)): u2 = Atn(Sqr(1 - e2) * Tan(B2)) dL = L2 - L1

sa1 = Sin(u1) * Sin(u2): sa2 = Cos(u1) * Cos(u2)