大地測量學(xué)課程設(shè)計--白塞爾大地主題解算

1、<p><b>　　大地測量學(xué)課程設(shè)計</b></p><p>　　設(shè)計題目：白塞爾大地主題解算</p><p>　　學(xué) 院： 礦業(yè)學(xué)院 </p><p>　　專 業(yè)： 測繪工程 </p><p>　　班 級： ******** </p>&l

2、t;p>　　學(xué) 號： * </p><p><b>　　目錄</b></p><p>　　1.基本原理及思想 1</p><p>　　2.白塞爾法大地主題正算步驟 2&

3、lt;/p><p>　　3.白塞爾法大地主題反算步驟 4</p><p>　　4.同一平行圈弧長、子午線弧長與大地線比較大小 6</p><p>　　5.程序代碼 8</p>

4、<p>　　6.演算示例 13</p><p>　　7.參考文獻(xiàn) 16</p><p>　　8.心得體會

5、17</p><p>　　9.教師評語 18</p><p><b>　　白塞爾大地主題解算</b></p><p><b>　　一：基本原理</b></p><p>　　建立以橢球中心為中心，以任意長（或單

6、位長）為半徑的輔助球，按以下三個步驟計算。</p><p>　　第一，   按一定條件將橢球面元素投影到輔助球面上。</p><p>　　第二，   在球面上解算大地問題。</p><p>　　第三，   將求得的球面元素按投影關(guān)系換算到相應(yīng)的橢球元素。</p><p>　　關(guān)鍵：確

7、定球面元素與橢球面元素的關(guān)系，即它們間的投影關(guān)系。</p><p>　　二：白塞爾法解算大地主題的基本思想:</p><p>　　以輔助球面為基礎(chǔ),將橢球面三角形轉(zhuǎn)換為輔助球面的相應(yīng)三角形,由三角形對應(yīng)元素關(guān)系,將橢球面上的大地元素按照白塞爾投影條件投影到輔助球面上，然后在球面上進(jìn)行大地主題解算，最后再將球面上的計算結(jié)果換算到橢球面上。</p><p>　　三：在球

8、面上進(jìn)行大地主題解算</p><p>　　球面上的大地主題正算：</p><p>　　球面上的大地主題反算：</p><p>　　球面三角元素間的相互關(guān)系：</p><p><b>　　四：正反算步驟</b></p><p>　　1.白塞爾法大地主題正算步驟</p><p>

9、;　　已知 、、、、()、S，計算、、()。</p><p>　?。?）將橢球面元素投影到球面上</p><p><b>　　由求：</b></p><p><b>　　計算輔助量和</b></p><p><b>　　,</b></p><p>　　計

10、算球面長度，將S化為</p><p><b>　　式中系數(shù)分別為：</b></p><p>　　上式右端含有，因此需要迭代計算。第一次迭代取近似值，第二次計算取</p><p>　　以后計算用代換代入上式迭代計算，直到所要求的精度為止。一般取。</p><p><b>　　解算球面三角形</b>&l

11、t;/p><p><b>　　計算</b></p><p><b>　　計算</b></p><p><b>　　或</b></p><p><b>　　計算</b></p><p>　?。?）將球面元素?fù)Q算到橢球面上</p>

12、;<p><b>　　由求</b></p><p><b>　　或</b></p><p>　　將球面經(jīng)差化為橢球面經(jīng)差l，求</p><p><b>　　l</b></p><p><b>　　式中</b></p><p

13、>　　式中的最大值為，故在計算時通?？梢月匀ゲ挥嫛?lt;/p><p><b>　　象限的判定</b></p><p><b>　　l</b></p><p><b>　　其中、為銳角。</b></p><p>　　2.白塞爾法大地主題反算步驟</p><

14、;p>　　已知 、、、、、，計算()、()、S。</p><p>　?。?）將橢球面元素投影到球面上</p><p><b>　　由B求u</b></p><p><b>　　，，l= </b></p><p><b>　　,</b></p><p&g

15、t;<b>　　,</b></p><p>　　采用逐次趨近法，由l計算</p><p>　　在反算中，已知橢球面上經(jīng)差l，球面經(jīng)差上的對應(yīng)經(jīng)差未知，為了由l求，由下式可知還需計算、、，計算又還需量，故需要進(jìn)行迭代計算。</p><p><b>　　第一次趨近，取l；</b></p><p><

16、;b>　　，</b></p><p><b>　　或</b></p><p><b>　　判斷的象限</b></p><p><b>　　判斷象限</b></p><p><b>　　l+</b></p><p>

17、　　仿照上述計算步驟迭代計算，直到為止。</p><p>　　將球面元素?fù)Q算到橢球面上</p><p><b>　　或</b></p><p>　　象限的判斷與前面一致</p><p>　　五：同一平行圈弧長、子午線弧長與大地線比較大小</p><p>　　子午線弧長計算公式：</p>

18、<p><b>　　式中：</b></p><p>　　平行圈弧長公式：（）</p><p>　　不同緯度對應(yīng)的一些弧長的數(shù)值</p><p>　　利用白塞爾大地主題反算求解大地線長S</p><p>　　緯度為30°，經(jīng)差為1°的平行圈弧長S=96 488m，兩點(diǎn)間大地線長為96 48

19、7.595m</p><p>　　經(jīng)度為30°，緯度差為1°的子午線弧長X=110 863m，兩點(diǎn)間大地線長為110 862.869m</p><p>　　通過比較可知，同一平行圈或同一子午線兩點(diǎn)間大地線長度與對應(yīng)的平行圈弧長或子午線弧長相等。</p><p><b>　　六：程序代碼</b></p><

20、p>　　#include<stdio.h></p><p>　　#include<math.h></p><p>　　double hudu(double,double,double); /*度分秒轉(zhuǎn)換為弧度*/</p><p>　　double du(double); /*弧度轉(zhuǎn)換為度*

21、/</p><p>　　double fen(double); /*弧度轉(zhuǎn)換為分*/</p><p>　　double miao(double); /*弧度轉(zhuǎn)換為秒*/</p><p>　　#define PI 3.1415926</p><p>　　void main

22、 (void)</p><p><b>　　{</b></p><p><b>　　int k;</b></p><p>　　printf("請選擇大地主題算法，若執(zhí)行正算，請輸入1；若執(zhí)行反算，請輸入2。\n");</p><p>　　scanf("%d",

23、&k);</p><p>　　/*大地主題正算*/</p><p><b>　　if(k==1)</b></p><p><b>　　{</b></p><p>　　double ax,ay,az,bx,by,bz,cx,cy,cz,S,dz,ez,fz,B1,B2,L1,L2,A1,A2;

24、</p><p>　　intdx,dy,ex,ey,fx,fy;</p><p>　　double e2,W1,sinu1,cosu1,sinA0,coto1,sin2o1,cos2o1,sin2o,cos2o,A,B,C,r,t,o0,o,g,sinu2,q;</p><p>　　/*輸入度分秒數(shù)據(jù)*/</p><p>　　printf(

25、"請輸入大地線起點(diǎn)緯度度分秒\n");</p><p>　　scanf("%lf%lf%lf",&ax,&ay,&az);</p><p>　　printf("請輸入大地線起點(diǎn)經(jīng)度度分秒\n");</p><p>　　scanf("%lf%lf%lf",&

26、bx,&by,&bz);</p><p>　　printf("請輸入大地方位角度分秒\n");</p><p>　　scanf("%lf%lf%lf",&cx,&cy,&cz);</p><p>　　printf("請輸入大地線長度\n");</p>&

27、lt;p>　　scanf("%lf",&S);</p><p><b>　　/*調(diào)用函數(shù)*/</b></p><p>　　B1=hudu(ax,ay,az);</p><p>　　L1=hudu(bx,by,bz);</p><p>　　A1=hudu(cx,cy,cz);</p&

28、gt;<p>　　/*白塞爾大地主題解算*/</p><p>　　e2=0.006693421622966;</p><p>　　W1=sqrt(1-e2*sin(B1)*sin(B1));</p><p>　　sinu1=sin(B1)*(sqrt(1-e2))/W1;</p><p>　　cosu1=cos(B1)/W1;&

29、lt;/p><p>　　sinA0=cosu1*sin(A1);</p><p>　　coto1=cosu1*cos(A1)/sinu1;</p><p>　　sin2o1=2*coto1/(coto1*coto1+1);</p><p>　　cos2o1=(coto1*coto1-1)/(coto1*coto1+1);</p>&

30、lt;p>　　A=6356863.020+(10718.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);</p><p>　　B=(5354.469-8.978*(1-sinA0*sinA0))*(1-sinA0*sinA0);</p><p>　　C=(2.238*(1-sinA0*sinA0))*(1-sinA0*sinA0)+0.006;

31、</p><p>　　r=691.46768-(0.58143-0.00144*(1-sinA0*sinA0))*(1-sinA0*sinA0);</p><p>　　t=(0.2907-0.0010*(1-sinA0*sinA0))*(1-sinA0*sinA0);</p><p>　　o0=(S-(B+C*cos2o1)*sin2o1)/A;</p>

32、<p>　　sin2o=sin2o1*cos(2*o0)+cos2o1*sin(2*o0);</p><p>　　cos2o=cos2o1*cos(2*o0)-sin2o1*sin(2*o0);</p><p>　　o=o0+(B+5*C*cos2o)*sin2o/A;</p><p>　　g=(r*o+t*(sin2o-sin2o1))*sinA0;

33、</p><p><b>　　/*求B2*/</b></p><p>　　sinu2=sinu1*cos(o)+cosu1*cos(A1)*sin(o);</p><p>　　B2=atan(sinu2/(sqrt(1-e2)*sqrt(1-sinu2*sinu2)));</p><p><b>　　/*求L2

34、*/</b></p><p>　　q=atan(sin(A1)*sin(o)/(cosu1*cos(o)-sinu1*sin(o)*cos(A1)));</p><p><b>　　/*判斷q*/</b></p><p>　　if(sin(A1)>0 && tan(q)>0)</p><

35、;p>　　q=fabs(q);</p><p>　　else if(sin(A1)>0 && tan(q)<0)</p><p>　　q=PI-fabs(q);</p><p>　　else if(sin(A1)<0 && tan(q)<0)</p><p>　　q=-fabs(q

36、);</p><p><b>　　else</b></p><p>　　q=fabs(q)-PI;</p><p>　　L2=L1+q-g/3600/180*PI;</p><p><b>　　/*求A2*/</b></p><p>　　A2=atan(cosu1*sin(A

37、1)/(cosu1*cos(o)*cos(A1)-sinu1*sin(o)));</p><p><b>　　/*判斷A2*/</b></p><p>　　if(sin(A1)<0 && tan(A2)>0)</p><p>　　A2=fabs(A2);</p><p>　　else if(s

38、in(A1)<0 && tan(A2)<0)</p><p>　　A2=PI-fabs(A2);</p><p>　　else if(sin(A1)>0 && tan(A2)>0)</p><p>　　A2=PI+fabs(A2);</p><p><b>　　else<

39、/b></p><p>　　A2=2*PI-fabs(A2);</p><p><b>　　/*調(diào)用函數(shù)*/</b></p><p>　　dx=(int)(du(B2));</p><p>　　dy=(int)(fen(B2));</p><p>　　dz=miao(B2);</p&g

40、t;<p>　　ex=(int)(du(L2));</p><p>　　ey=(int)(fen(L2));</p><p>　　ez=miao(L2);</p><p>　　fx=(int)(du(A2));</p><p>　　fy=(int)(fen(A2));</p><p>　　fz=miao(

41、A2);</p><p>　　printf("大地線終點(diǎn)緯度度分秒分別為：\n%d\n%d\n%lf\n",dx,dy,dz);</p><p>　　printf("大地線終點(diǎn)經(jīng)度度分秒分別為：\n%d\n%d\n%lf\n",ex,ey,ez);</p><p>　　printf("終點(diǎn)大地方位角度分秒分別為：\n

42、%d\n%d\n%lf\n",fx,fy,fz);</p><p><b>　　}</b></p><p>　　/*大地主題反算*/</p><p>　　else if(k==2)</p><p><b>　　{</b></p><p>　　double fx,fy

43、,fz,gx,gy,gz,hx,hy,hz,ix,iy,iz,jz,kz,B1,B2,L1,L2,S,A1,A2;</p><p>　　int jx,jy,kx,ky;</p><p>　　double e2,W1,W2,sinu1,sinu2,cosu1,cosu2,L,a1,a2,b1,b2,g,g2,g0=0.0,r,p,q,sino,coso,o,sinA0,x,t1,t2,A,B

44、,C,y;</p><p>　　/*輸入度分秒數(shù)據(jù)*/</p><p>　　printf("請輸入大地線起點(diǎn)緯度度分秒\n");</p><p>　　scanf("%lf%lf%lf",&fx,&fy,&fz);</p><p>　　printf("請輸入大地線起點(diǎn)經(jīng)度

45、度分秒\n");</p><p>　　scanf("%lf%lf%lf",&gx,&gy,&gz);</p><p>　　printf("請輸入大地線終點(diǎn)緯度度分秒\n");</p><p>　　scanf("%lf%lf%lf",&hx,&hy,&

46、hz);</p><p>　　printf("請輸入大地線終點(diǎn)經(jīng)度度分秒\n");</p><p>　　scanf("%lf%lf%lf",&ix,&iy,&iz);</p><p><b>　　/*調(diào)用函數(shù)*/</b></p><p>　　B1=hudu

47、(fx,fy,fz);</p><p>　　L1=hudu(gx,gy,gz);</p><p>　　B2=hudu(hx,hy,hz);</p><p>　　L2=hudu(ix,iy,iz);</p><p>　　/*白塞爾大地主題解算*/ </p><p>　　e2=0.006693421622966;</

48、p><p>　　W1=sqrt(1-e2*sin(B1)*sin(B1));</p><p>　　W2=sqrt(1-e2*sin(B2)*sin(B2));</p><p>　　sinu1=sin(B1)*sqrt(1-e2)/W1;</p><p>　　sinu2=sin(B2)*sqrt(1-e2)/W2;</p><p

49、>　　cosu1=cos(B1)/W1;</p><p>　　cosu2=cos(B2)/W2;</p><p><b>　　L=L2-L1;</b></p><p>　　a1=sinu1*sinu2;</p><p>　　a2=cosu1*cosu2;</p><p>　　b1=cosu1

50、*sinu2;</p><p>　　b2=sinu1*cosu2;</p><p>　　/*逐次趨近法求解A1*/</p><p><b>　　g=0;</b></p><p><b>　　r=L;</b></p><p><b>　　while(1)</b&

51、gt;</p><p><b>　　{</b></p><p>　　p=cosu2*sin(r);</p><p>　　q=b1-b2*cos(r);</p><p>　　A1=atan(p/q);</p><p><b>　　/*判斷A1*/</b></p>

52、<p>　　if(p>0 && q>0)</p><p>　　A1=fabs(A1);</p><p>　　else if(p>0 && q<0)</p><p>　　A1=PI-fabs(A1);</p><p>　　else if(p<0 && q&l

53、t;0)</p><p>　　A1=PI+fabs(A1);</p><p><b>　　else</b></p><p>　　A1=2*PI-fabs(A1);</p><p>　　sino=p*sin(A1)+q*cos(A1);</p><p>　　coso=a1+a2*cos(r);<

54、;/p><p>　　o=atan(sino/coso);</p><p><b>　　/*判斷o*/</b></p><p>　　if(coso>0)</p><p>　　o=fabs(o);</p><p><b>　　else</b></p><p&

55、gt;　　o=PI-fabs(o);</p><p>　　sinA0=cosu1*sin(A1);</p><p>　　x=2*a1-(1-sinA0*sinA0)*coso;</p><p>　　t1=(33523299-(28189-70*(1-sinA0*sinA0))*(1-sinA0*sinA0))*1e-10;</p><p>　

56、　t2=(28189-94*(1-sinA0*sinA0))*1e-10;</p><p>　　g2=(t1*o-t2*x*sino)*sinA0;</p><p>　　if(fabs(g2-g0)<=1e-100)</p><p><b>　　break;</b></p><p><b>　　else

57、</b></p><p><b>　　{</b></p><p><b>　　r=L+g2;</b></p><p><b>　　g0=g2;</b></p><p><b>　　}</b></p><p><b&

58、gt;　　}</b></p><p><b>　　/*求解S*/</b></p><p>　　A=6356863.020+(10708.949-13.474*(1-sinA0*sinA0))*(1-sinA0*sinA0);</p><p>　　B=10708.938-17.956*(1-sinA0*sinA0);</p>

59、;<p><b>　　C=4.487;</b></p><p>　　y=((1-sinA0*sinA0)*(1-sinA0*sinA0)-2*x*x)*coso;</p><p>　　S=A*o+(B*x+C*y)*sino;</p><p><b>　　/*求解A2*/</b></p><

60、;p>　　A2=atan(cosu1*sin(r)/(b1*cos(r)-b2));</p><p><b>　　/*判斷A2*/</b></p><p>　　if(p<0 && q<0)</p><p>　　A2=fabs(A2);</p><p>　　else if(p<0 &

61、amp;& q>0)</p><p>　　A2=PI-fabs(A2);</p><p>　　else if(p>0 && q>0)</p><p>　　A2=PI+fabs(A2);</p><p><b>　　else</b></p><p>　　A2

62、=2*PI-fabs(A2);</p><p><b>　　/*調(diào)用函數(shù)*/</b></p><p>　　jx=(int)(du(A1));</p><p>　　jy=(int)(fen(A1));</p><p>　　jz=miao(A1);</p><p>　　kx=(int)(du(A2))

63、;</p><p>　　ky=(int)(fen(A2));</p><p>　　kz=miao(A2);</p><p>　　printf("起點(diǎn)大地方位角度分秒分別為：\n%d\n%d\n%lf\n",jx,jy,jz);</p><p>　　printf("終點(diǎn)大地方位角度分秒分別為：\n%d\n%d\n%

64、lf\n",kx,ky,kz);</p><p>　　printf("大地線長度為：%lf\n",S);</p><p><b>　　}</b></p><p><b>　　/*數(shù)據(jù)錯誤*/</b></p><p><b>　　else</b>&l

65、t;/p><p>　　printf("數(shù)據(jù)錯誤，請重新輸入\n");</p><p><b>　　}</b></p><p>　　/*度分秒轉(zhuǎn)換為弧度*/</p><p>　　double hudu(double x,double y,double z) </p><p>&

66、lt;b>　　{</b></p><p>　　double A0;</p><p>　　A0=(x+y/60+z/3600)*PI/180;</p><p>　　return A0;</p><p><b>　　}</b></p><p>　　/*弧度轉(zhuǎn)換為度*/++++++++

67、++++++</p><p>　　double du(double B0)</p><p><b>　　{</b></p><p>　　double x0;</p><p>　　x0=(int)(B0*180/PI);</p><p>　　return x0;</p><p&g

68、t;<b>　　}</b></p><p>　　/*弧度轉(zhuǎn)換為分*/</p><p>　　double fen(double C0)</p><p><b>　　{</b></p><p>　　double _y,y0;</p><p>　　_y=(int)(C0*180/P

69、I);</p><p>　　y0=(fabs)((int)((C0*180/PI-_y)*60));</p><p>　　return y0;</p><p><b>　　}</b></p><p>　　/*弧度轉(zhuǎn)換為秒*/</p><p>　　double miao(double D0)<

70、/p><p><b>　　{</b></p><p>　　double _z1,_z2,z0;</p><p>　　_z1=(int)(D0*180/PI);</p><p>　　_z2=(int)((D0*180/PI-_z1)*60);</p><p>　　z0=(fabs)((double)((

71、(D0*180/PI-_z1)*60-_z2)*60));</p><p>　　return z0;</p><p><b>　　}</b></p><p><b>　　七：演算示例</b></p><p><b>　　正算</b></p><p>&l

72、t;b>　　反算</b></p><p><b>　　正算計算結(jié)果</b></p><p><b>　　反算計算結(jié)果</b></p><p><b>　　八：參考文獻(xiàn)</b></p><p>　　[1]史國友,周曉明,賈傳熒. 貝塞爾大地主題正解的改進(jìn)算法[J

73、]. 大連海事大學(xué)學(xué)報,2008,01:15-19. </p><p>　　[2]郭際明,丁士俊,蘇新洲,劉宗泉.大地測量學(xué)基礎(chǔ)實踐教程.武漢:武漢大學(xué)出版社,2009</p><p>　　[3]裴連磊. 用C語言實現(xiàn)大地主題解算[J]. 價值工程,2013,20:235-236. </p><p>　　[4]周振宇,郭廣禮,賈新果. 大地主題解算方法綜述[J].

74、測繪科學(xué),2007,04:190-191+174+200. </p><p>　　[5]丁士俊,楊艷梅,史俊波,程新明. 大地主題解算幾種不同算法在計算中應(yīng)注意的問題[J]. 黑龍江工程學(xué)院學(xué)報(自然科學(xué)版),2013,03:1-5. </p><p>　　[6]王建強(qiáng),胡明慶. 貝賽爾大地主題解算分析[J]. 測繪科學(xué),2012,01:30-31. </p><p&g

75、t;　　[7]徐曉晗,謝云開,李亞軍. 大地主題解算實用算法[J]. 科學(xué)技術(shù)與工程,2012,09:2062-2068. </p><p>　　[8]許厚澤. 關(guān)于正反大地主題解算方法的綜合研究[J]. 測量制圖學(xué)報,1958,04:274-288. </p><p><b>　　. </b></p><p><b>　　心得體會&

76、lt;/b></p><p>　　學(xué)習(xí)課程設(shè)計是培養(yǎng)學(xué)生綜合運(yùn)用所學(xué)知識,發(fā)現(xiàn),提出,分析和解決實際問題,鍛煉實踐能力的重要環(huán)節(jié),是對學(xué)生實際工作能力的具體訓(xùn)練和考察過程.隨著科學(xué)技術(shù)發(fā)展的日新日異，測繪技術(shù)也不斷發(fā)展，這次課程設(shè)計讓我進(jìn)一步了解和掌握大地線是橢球面上最短程曲線的特性以及大地主題解算在測繪學(xué)科中的地位和在實際應(yīng)用中的意義，同時，我也學(xué)會利用計算機(jī)解決問題的基本技能得到訓(xùn)練和提高，為我走向工作

77、崗位打下一定實踐基礎(chǔ)。</p><p>　　這次課程設(shè)計最難得地方就是程序設(shè)計，雖然之前學(xué)過C語言，但是只是懂得一些基礎(chǔ)，而且之前學(xué)習(xí)了的c語言也忘記的差不多了，所以我又重新開始看書，復(fù)習(xí)了一下c語言才開始做課程設(shè)計在這次課程設(shè)計中，我發(fā)現(xiàn)教材上有一些印刷錯誤，雖然問題不大，但對于還不熟悉這門課的我們來說，在學(xué)習(xí)過程中難免會造成疑惑，多走彎路。此外，教材上的算法雖然簡單易懂，但是它在計算系數(shù)時是用橢球參數(shù)直接代入

78、后化簡的公式，一個程序只能依據(jù)一個參考橢球，不夠靈活。我在此基礎(chǔ)上稍微做了修改，可以以不同參考橢球為依據(jù)，使程序更具普遍性，讓課程設(shè)計更完善。</p><p>　　回想起起此次課程設(shè)計，至今我仍感慨頗多，從理論到實踐，在這段日子里，可以說得是苦多于甜，但是可以學(xué)到很多很多的東西，同時不僅可以鞏固了以前所學(xué)過的知識，而且學(xué)到了很多在書本上所沒有學(xué)到過的知識。通過這次課程設(shè)計使我懂得了理論與實際相結(jié)合是很重要的，只有

79、理論知識是遠(yuǎn)遠(yuǎn)不夠的，只有把所學(xué)的理論知識與實踐相結(jié)合起來，從理論中得出結(jié)論，才能真正學(xué)以致用，從而提高自己的實際動手能力和獨(dú)立思考的能力。在設(shè)計的過程中遇到問題，可以說得是困難重重，但可喜的是最終都得到了解決。</p><p>　　雖然課程設(shè)計已經(jīng)結(jié)束了，但是在這次的課程設(shè)計中不僅檢驗了我所學(xué)習(xí)的知識，也培養(yǎng)了我如何去把握一件事情，如何去做一件事情，又如何完成一件事情。學(xué)會了合作，學(xué)會了運(yùn)籌帷幄。課程設(shè)計是我們