椭球大地测量学

1、614贝塞尔微分方程的其他解法推导贝塞尔大地问题解算公式，关键在于求解贝塞尔微分方程：dS 二。*1 e2cos2 ud ，dl =- e2cos2 udk（1）解算这组微分方程要用到椭圆积分。而椭圆积分的原函数一级不能用初等函数表示；贝塞尔 将被积函数展开，逐项积分，把原函数表示为具有一定精度的比较简单的表达式。此外，还 有其他各种不同的解法：或者将微分方程的形式加以改变，或者采用不同的积分方法。因而 出现了各式各样的公式。详细评述各种解法是因难的。下面只对一些主要公式作简单的介 绍。了解前人的工作，可能对后来的工作有所启发。一、维罗维茨方法利用d。、d入与du的关系，将（1）式改化为 TO

2、C o 1-5 h z : cl 一 e2 cos2 u cosuc%1 e2 cos2 u .dS =:du， dl =, du（1）cos2 u 一 c2cosu、：cos2 u 一 c2式中c = cos u - sin A（2 ）然后将被积函数展开，逐项积分。公式结构比较复杂55。二、安德列也夫方法将（1）式改化为IdS = a1：忐 dA， dl =如气1-一L dAsin Asin2 A c2vsin2 A c2同样展开被积函数积分，公式也很复杂【278。三、鼓尔默持法赫尔默特同样采用贝塞尔辅助球，但规定大地线及 其投影后的大圆弧从其纬度最高点起算，向东为正；如 图 610。由图知

3、sin u = sin u cos b0cos2 u = 1 sin2 u （1 sin2 b）代入（1）第一式，得e2 sin2 uds = a. 1e2 cosu ,1sin2 bdb01e2 cosu10引入辅助量K12:K2 = e2 sin2 B0代入前式，经变换，可得dS = a e2 J + K; + 2K;cos2bdb1同样对于（1）式的第二个微分方程有dl = d人一 cosu (1 - 1 - e2)+ K12cos2。d。01 一 sin 2 u cos 2。然后将上式展开积分。赫尔驮特采用级数回求法，在由s求。时，消除了迭代计算。用予正解比较方便，而反 解仍需迭代。1

4、954年博德米勒曾依据海福特椭球将贝塞尔赫尔默特公式各系数编成了算表。1957 年史腊德又对该算表作了改进。1974年总参测绘局第二测绘大队编制过依据克拉索夫基斯 椭球的算表。四、韦贝尔方法韦贝尔在贝塞尔和赫尔默特方法的基础上，作了一些改进。他把子午线弧长的计算公式 转化为大地线的长度计算，并且采用具有一定精度的封闭公式，结构比较简单。同赫尔默特 方法一样，在正解中由s求。无需达代，而反解仍需迭代。详细说明可参考13。五、勒瓦路易一一杜皮方法参考 6-8，有 TOC o 1-5 h z s 1.，、1.，dS = b1 + K2 sin2（M +。）一 K4 sin4（M +。）+d。28-j

5、 sin( M +b) d。= J84j sin2（M +。）d。= J，积分得S = b。+ K2AJ2 - K4AJ4 + （7）同样可以改化经差计算公式。J2、J4人可由瓦利斯积分数值表中查取。瓦利斯积分表是法国人编算的，间隔太大，内 插不甚方便，计算精度不太高。因此，这种方法没有在大地测量中得到广泛的应用。六、约尔旦方法已知V = 1/、1 一 e2 cos2 u贝塞尔微分方程（1）式可以写成ad。= VdS，dk = Vdl（1）约尔旦方法仍以贝塞尔微分方程为基础，但是，不是将v展为e2或e2的幕级数逐项积 分，而是利用台劳级数值直接将。展为AS的幕级数，将入展为】的幕级数。它们的一

6、般 形式是：(8)式中下标“1”表示这些倒数应依据P1点的B1和A1计算。由（1）式知道d。a = VdSd 2。dV dV dBa=dS 2dS dB dSdld 2 k _ dV _ dV dBdl 2 dl dB dl(9)注意到dV 门 2dB V 3=t， = cos A dB VdS c可得dVV2=-n 2 t cos A dSc进而有d 2VV 41一 dS = F 2 乜2 A(1 一 12 +n 2 3n 21 ) 一 12 sin2 A又因dB dB dS “=V 2 ctgA cos Bdl dS dl得到dV=-n 2Vctg sin B dl进而d 2Vcos2 B

7、 I1=一门2VCos2 A(1 -3t2 +门2 3q212) 一 12 sin2 Adl2sin 2 A将求得的各阶导数代入（9）式，再代入（8）式，即可实现。与S、人与l的变换。显然，约尔旦方法在使用上受到距离的限制。为了便于实际应用，约尔旦公式还要作一 些技术性处理，见51。七、数值积分法解算贝塞尔微分方程在于求解积分 TOC o 1-5 h z S =。气1一e2 cos2 ud。 ，l = j1-e2 cos2 udX（10）现将第二式加以改变：1= j。1 - （1 - % 1 - e 2 cos 2 u） d人e2 cos2 u0= X-fa3d人0 1 + 与； 1 e2 c

8、os2 u由 6-8（19）式知sin m 7d 人=d。cos2 u代入前式，得l =人一e2 sinmj。】d。（11）数学文献上讲述许多数值积分方法。文献290认为高斯数值积分法所需的节点最少， 精度最高。依据高斯数值积分公式b f (x)dx =(b - a) R f (x ) i i ai=1x = a + (b - a)V式中n是节点数，X是第i个节点的自变量，f( xi)是第i个节点处的函数值。R,和V是高 斯积分公式的常数，见附录七。1将高斯积分公式用于(10)和(11)式Jb V1 -e2 cos2 udb =b0Jb 1 db =b Ri =bEo 1 + J1 - e2

9、cos2 u .=1 1+ J1 - e2 cos2 u 1得 R X：1 - e 2 cos2 u = bEi=1式中于是有(13)=R ：1 - e2 cos2 ui=1=(14)i=1I1 + J1 - e2 cos2 uS = ab Ebl = X - e 2b sin mEX(15)X = l + e 2b sin mE人或者写成(16)baEb为了计算E。、E人，需要计算出节点Pi处的。由图6 - 11可以写出：sin u = sin u cos(V b) + cos u sin(V b )cos Ai1i1 ii(17 ) 据文献731研究，如果取五个节点，对于任意长 距离，长度

10、计算误差为0.001米，经纬度、方位 角的计算误差为0.0001 。八、嵌套系数法68导出的(10)式、(12)式和(20)式中，含有A、B、C、a、8、丫以及a 、8 等系数。而这些系数中包含有K2 ( =e2cos2m )和 K2( = e2cos2m)的各次幕。为了便于在电子计算机上计算，在保持原有精度的条件下(仍取 至e4项)，将公式加以调整。由 68(10)式S = bAb - Bsinb cos(2M +b) - Csin 2b cos(4M + 2b) TOC o 1-5 h z .一 131,一3、A = 1 + K2 -K4 = K2(1 - K2)464416111 八 1

11、、B = K 2-K4 = K 2(1 - K 2)41644y=C =上 K 4 =PiA 1288前式可以写成S = bAb - p sin b cos(2M + b)-史 sin b cosb cos(4M + 2b)4=bAjb p sinbcos(2M +b) -4cosb cos(4M + 2b)(18)同样，对于经差变换公式，由68 (20)式有l =人一 sin mlab + p sinb cos(2M +b)a=竺 + 竺竺 cos2 m，816，p=竺cos2 m = *cos2 m164可以写成以l = Xasinm b + 一cos2msinbcos(2M +b)(19

12、)(1)依据(1)式，有cos B = cos B sin a = cos B sin a(2)上述类型的公式有的文献叫嵌套系数公式文献262和136给出扩展至e 6项嵌套公式，并在其他细节上给予不同的说明。 嵌套系数公式与68给出的相应公式精度相同，但结构紧凑些，适于电算。当然，应用恢套公式仍然需要迭代。615保持纬度不变的大地投影一张志新公式贝塞尔投影是保持方位角不变的大地投影。张志新在推导长距离大地问题解公式时， 采用了保持纬度不变的大地投影。椭球面上的大地线投影到辅助球面上仍为大圆弧，但球面 纬度和方位角度按下式确定183184】185：sin a = cos B/cosBBn是椭球面

13、大地线最高点的大地纬度(90 -m)，90 -Bn相当于贝塞尔公式中的m。n (1)式给出了球面上的三个元素，还需要推求另外三个元素a 2、。和入。首先讨论a 2和A2的转换关系 将克莱劳定理用于球面上的大圆弧：cos 甲 =cos 甲 sin a = cos 甲 sin a(3)sin a = cos B / cos B另外，由循球面上大地线的克架劳定理: TOC o 1-5 h z cos u = cos u sin A = cos u sin A(4)n1121(4)式除(3)式，得cos Bcos B sin aC12)(5)cos ucos u sin A顾及cos B = V 1

14、e 2 cos u得至。sin a = sin A1或者顾及(4)式可得ctga = V ctgA , ctga = V ctgA(6)张志新投影方法同贝塞尔投影方法相比，纬度计算简单了，但方位角的计算复杂了。 下面讨论ds和d。、以与dA的微分方程。球面上大圆的微分方程为d里=cos ada，依据(1)式，上式可写成dB = cos ada，而椭球面上大地线微分方程为dB = cos AdS ， M比较以上两组方程，可以得到dk = sin a sec Bda小1，dL = sin A sec BdS Ncos a 7dS = Mda，cos AdL =竺些dkN ctgA(7)上式就是边长

15、和经差的微分方程。类似于贝塞尔的方法，对上式进行积分，即可求得S与、l与人的关系式:a = % S + P sin a cos(2M +a) + P Lin2M sin 2 M sin 2(M +a )sin2(M +a)(8) TOC o 1-5 h z W 23nk = l + y a sin my sin m sin a cos(2M +a)式中PP =1 b(3P = p(3913)1 一一e2 cos2 m e4 cos4 m e6 cos6 m464256)i913)e2 cos2 m e4 cos4 m e6 cos6 m464256)(9)P =p卫e4 cos4 m + 兰e

16、6 cos6 m36496)1353y = e 2 e 4(2 cos2 m) + e 6(1 cos2 m + cos4 m) 1 2161683,5y = p cos2 m上式展开至e6项，同贝塞尔公式展开至e6项的精度相当。616椭球面对球面的正形投影椭球面对球面的正形投影是数学制图学中的一个基本问题。这里用它来解算大地问题。 椭球面正形投影到球面上的方法，是高斯1844年提出的。下面说明这种投影的要点。 设辅助球的半径为R，球面经差为入，球面纬度为轩可以写出于午圈投影长度比Ml = Rd? /MdB（1）平行圈投影长度比mB = R cos pd! / N cos Bdl（2）依据正形

17、投影的理论（详细说明见第十章），投影长度比与方向无关，即 m为任意方向投影长度比。将上式代入（1）式、（2）式，得RdpR cos p d!MdBN cos B 矛（3）式可以写成工 dp = -dB 些 R cos pN cos B dldq q= R dpR cos pdq =dBN cos Bdq = dq ” dl（3）（4）（5）（6）d d! = Al进而由（6）式，得dq = Adq（9）现提出以下条件；椭球面子午线投影到球面仍为子午线，即球面经度仅为大地经度的函 数；椭球面平行圈投影到球面仍为平行圈，即球面纬度仅为大地纬度的函数。因此，可设（7） （8）（10）q = Aq +

18、 KA、K是投影常数。（8）式和（10）式建立起椭球面上点位（q、l）与球面上相应点位（q、人）的投影关系。q 叫等量纬度。至于等量纬度与对应的纬度之间的关系，可按以下方法确定：(1 - e 2)dB_jB MdB _jB (1 - e 2)dB _jB (1 - e 2 sin2 B e 2 cos2 B)dB O N cos B O (1 e 2 sin 2 B)cos B O (1 e 2 sin 2 B) cos BJ B dB q JO cos Be 2 cos BdB1 e 2 sin 2 B为便于计算第二个积分，设两端微分:e sin B = sin 中e cos BdB = c

19、os 中卯因此有q = J b 笠-e J5O cos B O cos q(B q = ln tg 45。+ 一k2-ln tge 45。+皂k 2tge 45。+ k 2=tge1 - cos(90 + q)1 + cos(90 + q)(1 + e sin B 丫 2=k 1 - e sin B )因此q = ln tg 45。+kB Y1 - e siiB 丫2 2 人 1 + e siiB)(11)(11)式是大地纬度B与椭球面等量纬度q之间的关系式。对于球面，因为e = 0,球面纬度与球面等量纬度q的关系式为(12)/ = ln tg 45+ k 2)在上述投影过程中，除了辅助球半径

20、R以外，还有两个待定常数A和K。确定这三个 常数，需要三个条件。前面讲到，正形投影的长度比与方向无关，仅与点位(纬度)有关：(13)m = m(B) = mB + (B - B )i展开上式:1 ,、+ _(B - B2(d2 m，-k dB 2)11 z 、+ ,(B - B )3d3 m-I dB 3 ) i(14)现提出一下三个条件：+为了使投影长度比接近于1，且随纬差的增大而作缓慢的变化，(15)(16)(17)(d2 m )-k dB 2 )1上式表明，在纬度为B的平行圈上，长度比等于1。这个平行圈叫做标准平行圈。从 标准平行圈向南向北延神，长度比的变化十分缓慢，仅与m对B的三阶导数

21、有关。将(7)式代入(3)式，得由此得dmdB顾及（16）式，有Rd甲R cos甲m aMdBN cos B卯_ a M cos甲dNcodB,Jsin甲d甲M sin B cos甲)AR +-NcosBdBN2 cos2 B J,J , M cos甲 sin 甲 M sin B cos甲)AR A *+VN2 cos2 BN2 cos2 B(19)(20)A M cos 甲 sin 甲N2 cos2 B+ M cos B sin 甲 N 2 cos2 B(2 1)+ cos B即A sin 中sin B为了满足(17 )式，可以要求(- A sin 中 + sin B) 0 dB而dd中M

22、cos2 中(A sin q + sin B) A cos q + cos B A2dBdBN cos B因而有A2 M1 cos2 q1- cosB(22)N cos B1或A22 N cos2 B(1 e2 sin2 B )cos2 BM1 e 2由(21)式有A 2 A 2 sin 2 q _ A 2 sin 2 B因此A2. 2 b +(1 e2 sin2 B )cos2 B】+ e2 cos4 B11 e 21 e 2于是确定了常数 A：A ：1 + e2 cos4 B（23）常数K按下述方法确定，由B和A按（21）式计算。依据B由（11）式、（12）式计算q1和q；。按（10）式求

23、K：(24)最后，由条件刀广1，写出A R cosR = 1N1 cos B顾及（22）式，有R = N cos B / A cos 中=N cos B /N / M cos B = J：M N（25）i ii i i i i ii i由此得出，辅助球的半径等于标准纬度处的平均曲率半径。至此三个待定常数已经完全确定了。于是通过（7）式可由l求得入，通过（i0）式可由 q求得q。而由（ii）式、（i2）式知道，q是B的函数，q是p的函数。因此，椭球面 上的一定点（B、1），在球面上有唯一的点（p、”与之相应。下面讨论长度比的变化。略去推导，约出长度比对纬度的三阶导数：d 2 m2e 2（i e

24、2 ）sin 2 B=-q 2e 2 sin 2BdB 2（i e 2 sin2 B）2代入（i8）式，得m = i msin2B（B B2（26）取 B=45，B- B=i.5，则m i = sinZBB B2 q i -10-8由此知道，当纬差小于3度时，长度变形（mi）的最大值小于ii0-8。这个误差通 常可以不计。因此，对于纬幅小于3度的地区，如果以平均纬度处的平均曲率半径为半径作 辅助球，那么，可用辅助球面代替椭球面。两个面上对应的角度相等，长度变形可以忽略。就是说，在一定的范围（AB - 3 ）和精度要求（m i - i -i0 8）之内，可将椭球面上的 大地问题解算转化到球面上来

25、解算。6-17应用椭球面对球面的正形投影解算大地问题（一）高斯平均引数公式的另一推导方法槌球而球而99 6-13村球面正推毁曜到球3L现在利用上节的结果解算大地问题。1.椭球面元素正形投影到球面上设椭球面上PP2两点的大地纬度为BB2,大地经度为L、L2,两点间的距离为S, 前进大地方位角为AAJ。经差、纬差、方位角差为(1)(3)l = L - L , b = B - B , a = A - A TOC o 1-5 h z 212121又设中纬度、中方位角为B = !(B + B )(2), A = !(A + A)m 212m 212将椭球面三角形P.PN正形投影到半径为R =，、jM N

26、的辅助球面上，得PPN。1 2m m m1 2P、P2的球面纬度为、,经差为人。设甲 =2G +甲)，甲=甲一甲(4)由上节(8)式知道人=Al(5)由上节(22)式知道A =二土(6)cos 中 R而PP的方位角与PR的方位角相等。因为长度变形可以忽略，因此，PP=b，而 1 21 21 22.解算球面三角形利用球面三角的高斯一一德兰布尔公式，顾及(1)、(2)、(3)、(4)各式，可以写出: 5 . 4人，cQ力Asin sin A=sin cos ffi,sin cos A=cos sin12m2m2m22Q a人小q a力A甲卜cos sin =sin sin ffi,cos cos

27、=cos sin(7)222m2222 TOC o 1-5 h z 将上式中的小角度正弦、余弦函数展开，得J人2 )sin Am=X 1 -一 cos 中24 )cos Am丫1- /(8)a21 24八r1VY1 -7V87最后一式可以写为J(9)由(8)式前三式略去高次项，化简得到X = b sin A sec 9中=b cos Ama = b sin A tg 中 1 +(b 2X 2 k1 药 F J(i X2a 2 )1 + 一+ 一 HYPERLINK l bookmark31 o Current Document k1224 J/ b 2 a 2 )+ 一24 /(10)123.

28、球面元素转化到椭球面上将（5）式、（6）式代入（10）式的第一式另外，由于代入（10）式的第二式，得 sin Asec B 1 一 + k 24 24)mm人Sb =Mm(12 cos A 1 +a2+12 24/(11)(12)(13)最后，方位角公式可以写为Sa =Nmsin A tgB b 2 a 21 + T2 + 24 J(14)顾及b 2 a 2b 2 + a 2 a 2+ =12 2412_ b212 a2 _ b2 a212 cos2 B=+=+m24 12 12 24 12 2412最后得到S Al =sin AN mmsec B k1+命 - /人S Ab =cos AM

29、mmC a 2 12)1 + + k 24 12J(15)S , a2 , 12 cos2 B , b2）,a =sin A tgB I 1 + +m + /Nmm m 241212 上式与 6-4 （ 17）式相同。这里应用完全不同于6-4的方法，又一次导出了高斯平均引数大地问题解算公式。可 以看到，在高斯解法中，不是直接特大地纬度化为球面纬度，而是运用根简单的方法将球面 纬差化为椭球面纬差，避开了等量纬度，计算大为简易。618应用椭球面对球面的正形投影解算大问题（二一博林公式1981年博林应用椭球面对面的正形投影，导出了一组大地问题解算公式202。与高斯解 法不同（高斯取平均纬圈为标准平行

30、圈），博林取通过P1点的平行圈作为标准平行圈，因此 在正解中无需迭代。对于短距离（SV150公里）大地问题，仍能保持足够的精度。博林公 式的推导有某些特点，结构不算复杂。文献232曾推荐博林公式作为短距离大地问题解算的实用公式。下面依据随林的基本思想，采用本书惯用的错号和已有公式，详细叙述解算方法。公式 推演过程有所简化，而最终形式基本与原著相同。一、大地问题正解仍然使用图6-13，有以下关系吼=cos p sin Acosp cosb cos A - sin p sin b（1）（2）S ，一，.(1 + e 2 cos2 B )R1。圮1 + e 2顾及sin 中=-sin B1 A 1（

31、6-16（21）式）（3）cos p =、；1 - sin 2 p =.1 -sin2 B =tcosB: A 21 A 1（4）A = .,-1 + e2 cos4 B 1(3)、(4)式代入(1)式，得式中I匕=：1 + e，2 cos 2 B另外有（3）式、（4）sin A1AV cos B2 V cos B cos A cos b - sin B sin bi i ii匕 sin Acosb (V cos A - tgB tg。)tgX =c o 甲 c o b - s i np s i nr c o A式代入上式，顾及人=Al得，tgAl =Atgb sin A1V cos B -

32、sin B tgb cos A1111Atgb s i nAV co B - s i nB tgb co A1 -1（5）（6）（7）为了计算P2、与之间的纬差，利用球面子午弧长AX，等于椭球面对应的子午线弧长AX的关系，有X =X = R 中（8）式中平=甲2-p1由3-7（12）式知e 2cos82 B AB 2m7B)cos2B B 2m式中NB = B2 - B，AB B = B +121依据附录八展开Mm3n211 -1-V 21NB代入（9）式，并且注意1/M1=VJ %略去高次项，得AD N XN B =R1cos2B NB2 1 -NB - -n 2(1 - 12)NB22 V 28i= NXvR 11顾及（8）式，经过整理，、cos 2B NB28 m3 e2y NB - 3 e, 2 cos 2 B NB 2V281保持原有精度，NB = Np V1 1 -e 2cos 2B NB 221 NB V2 e 2 sin 2B 因为A “ 13 e2NB (2=Np