攝影測量實驗報告(空間后方交會—前方交會)(共11頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上空間后方交會-空間前方交會程序編程實驗一 實驗?zāi)康囊笳莆者\用空間后方交會-空間前方交會求解地面點的空間位置。學(xué)會運用空間后方交會的原理，根據(jù)所給控制點的地面攝影測量坐標(biāo)系坐標(biāo)以及相應(yīng)的像平面坐標(biāo)系中的坐標(biāo)，利用計算機編程語言實現(xiàn)空間后方交會的過程，完成所給像對中兩張像片各自的外方位元素的求解。然后根據(jù)空間后方交會所得的兩張像片的內(nèi)外方位元素，利用同名像點在左右像片上的坐標(biāo)，求解其對應(yīng)的地面點在攝影測量坐標(biāo)系中的坐標(biāo)，并完成精度評定過程，利用計算機編程語言實現(xiàn)此過程。二 儀器用具 計算機、編程軟件（MATLAB）三 實驗數(shù)據(jù)實驗數(shù)據(jù)包含四個地面控制點（GCP）的地面攝

2、影測量坐標(biāo)及在左右像片中的像平面坐標(biāo)。此四對坐標(biāo)運用最小二乘法求解左右像片的外方位元素，即完成了空間后方的過程。另外還給出了5對地面點在左右像片中的像平面坐標(biāo)和左右像片的內(nèi)方位元素。實驗數(shù)據(jù)如下：點號 左片 右片 地面攝影測量坐標(biāo) x y x y X Y Z GCP1 16.012 79.963 -73.93 78.706 5083.205 5852.099 527.925 GCP2 88.56 81.134 -5.252 78.184 5780.02 5906.365 571.549 GCP3 13.362 -79.37 -79.122 -78.879 5210.879 4258.446 4

3、61.81 GCP4 82.24 -80.027 -9.887 -80.089 5909.264 4314.283 455.484 1 51.758 81.555 -39.953 78.463 2 14.618 -0.231 -76.016 0.036 3 49.88 -0.792 -42.201 -1.022 4 86.243 -1.346 -7.706 -2.112 5 48.135 -79.962 -44.438 -79.736 內(nèi)方位元素：f=152.000mm，x0=0，y0=0四 實驗框圖輸入GCP的像點坐標(biāo) xy確定初始值=0，Xs，Ys，Zs計算旋轉(zhuǎn)矩陣R計算像點在像空間坐標(biāo)系

4、中的近似值（x）,（y）,并組成誤差方程的常數(shù)項l計算誤差方程的系數(shù)項組成系數(shù)矩陣A組成法方程式，計算系數(shù)AA常數(shù)項AL求解外方位元素計算、Xs、Ys、Zs改正后的值改正數(shù)是否小于限差？ 大于小于計算完畢此過程完成空間后方交會求解像片的外方位元素，其中改正數(shù)小于限差（0.00003，相當(dāng)于0.1的角度值）為止。在這個過程中采用迭代的方法，是外方位元素逐漸收斂于理論值，每次迭代所得的改正數(shù)都應(yīng)加到上一次的初始值之中。輸入所需計算點的像平面坐標(biāo)x1,y1；x2,y2根據(jù)后方交會所得的旋轉(zhuǎn)矩陣R1,R2計算像點在左右像空間輔助坐標(biāo)系中的坐標(biāo)X1 Y1 Z1，X2 Y2 Z2 計算攝影基線的三個坐標(biāo)分

5、量Bx By Bz 計算個點在左右像片中的的投影系數(shù)N1 N2計算地面所求點在地面攝影測量坐標(biāo)系中的坐標(biāo)XA YA ZA計算完畢，精度評定在空間后方交會中運用的數(shù)學(xué)模型為共線方程確定Xs，Ys，Zs的初始值時，對于左片可取地面左邊兩個GCP的坐標(biāo)的平均值作為左片Xs 和Ys的初始值，取右邊兩個GCP的坐標(biāo)平均值作為右片Xs 和Ys的初始值。Zs可取地面所有GCP的Z坐標(biāo)的平均值再加上航高。空間前方交會的數(shù)學(xué)模型為： 五 實驗源代碼function Main_KJQHFJH()global R g1 g2 m G a c b1 b2;m=10000;a=5;c=4;feval(shuru); %

6、調(diào)用shuru（）shurujcp（）函數(shù)完成像點及feval(shurujcp); %CCP有關(guān)數(shù)據(jù)的輸入 XYZ=feval(MQZqianfangjh); %調(diào)用MQZqianfangjh（）函數(shù)完成空間前方、% 單位權(quán)中誤差 % %后方交會計算解得外方位元素global V1 V2; %由于以上三個函數(shù)定義在外部文件中故需VV=; %用feval（）完成調(diào)用過程for i=1:2*c VV(i)=V1(i);VV(2*i+1)=V2(i);endm0=sqrt(VV*(VV')/(2*c-6);disp('單位權(quán)中誤差m0為正負(fù):');disp(m0); %計算

7、單位權(quán)中誤差并將其輸出顯示輸入GCP像點坐標(biāo)及地面攝影測量坐標(biāo)系坐標(biāo)的函數(shù)和輸入所求點像點坐標(biāo)函數(shù)：function shurujcp()global c m;m=input('攝影比例尺：'); %輸入GCP像點坐標(biāo)數(shù)據(jù)函數(shù)并分別將其c=input('GCP的總數(shù)='); % 存入到不同的矩陣之中disp('GCP左片像框標(biāo)坐標(biāo)：');global g1;g1=zeros(c,2);i=1;while i<=c m=input('x='); n=input('y='); g1(i,1)=m;g1(i,2)=

8、n; i=i+1;end disp('GCP右片像框標(biāo)坐標(biāo)：');global g2;g2=zeros(c,2);i=1;while i<=c m=input('x='); n=input('y='); g2(i,1)=m;g2(i,2)=n; i=i+1; end%function shuru() global a;a=input('計算總像對點數(shù)='); %完成想計算所需的像平面坐標(biāo)global b1; %坐標(biāo)輸入，存入不同的矩陣中b1=zeros(a,2);disp('左片像點坐標(biāo)：')i=1;whil

9、e i<=a m=input('x='); n=input('y='); b1(i,1)=m;b1(i,2)=n; i=i+1;end %global b2;b2=zeros(a,2);disp('右片像點坐標(biāo)：')i=1;while i<=a m=input('x='); n=input('y='); b2(i,1)=m;b2(i,2)=n; i=i+1;end %global c;c=input('GCP的總數(shù)=');disp('GCP攝影測量系坐標(biāo)：')global

10、 G;G=zeros(3,c);i=1;while i<=c m=input('X='); n=input('Y='); v=input('Z='); G(i,1)=m;G(i,2)=n;G(i,3)=v; i=i+1;end%空間前方交會和后方交會函數(shù)：function XYZ=MQZqianfangjh()global R1 R2 a f b1 b2 Ra Rb;global X1 X2;R1=Ra;R2=Rb;R1=zeros(3,3);R2=zeros(3,3);global g1 g2 G V1 V2 V WF c QXX QXX

11、1 QXX2;xs0=(G(1,1)+G(3,1)/2;ys0=(G(1,2)+G(3,2)/2;Xs1,Ys1,Zs1,q1,w1,k1 R=houfangjh(g1,xs0,ys0); %對左片調(diào)用后方交會函數(shù)R1=R;V1=V;WF1=WF;QXX1=QXX;save '左片外方位元素為.txt' WF -ascii %將計算所得的外方位元素存入到.txt % 文件中for i=1:c g1(i,1)=g1(i,1)+V1(2*i-1);g1(i,2)=g1(i,2)+V1(2*i);endsave '左片像點坐標(biāo).txt' g1 -ascii xs0=(

12、G(2,1)+G(4,1)/2;ys0=(G(2,2)+G(4,2)/2;Xs2,Ys2,Zs2,q2,w2,k2 R=houfangjh(g2,xs0,ys0); %對右片調(diào)用后方交會函數(shù) R2=R; V2=V;WF2=WF;QXX2=QXX;save '右片外方位元素為.txt' WF ascii %將計算所得的外方位元素存入到.txt % 文件中for i=1:c g2(i,1)=g2(i,1)+V2(2*i-1);g2(i,2)=g2(i,2)+V2(2*i);endsave '右片像點坐標(biāo).txt' g2 -ascii X1=zeros(a,3);X2

13、=zeros(a,3);xx=zeros(3,1);xxx=zeros(3,1);for i=1:a ss=b1(i,1);b1(i,2);-f;dd=b2(i,1);b2(i,2);-f; xx=R1*ss; X1(i,:)=xx'xxx=R2*dd;X2(i,:)=xxx'endglobal Xs1 Xs2 Ys1 Ys2 Zs1 Zs2; BX=Xs2-Xs1;BY=Ys2-Ys1;BZ=Zs2-Zs1;global N1 N2;N1=zeros(1,a);N2=zeros(1,a);for i=1:a N1(1,i)=(BX*X2(i,3)-BZ*X2(i,1)/(X1

14、(i,1)*X2(i,3)-X2(i,1)*X1(i,3); N2(1,i)=(BX*X1(i,3)-BZ*X1(i,1)/(X1(i,1)*X2(i,3)-X2(i,1)*X1(i,3);end %計算投影系數(shù)，并計算五點的三維坐標(biāo)global XYZ;XYZ=zeros(a,3);for i=1:a XYZ(i,1)=Xs1+N1(1,i)*X1(i,1); XYZ(i,3)=Zs1+N1(1,i)*X1(i,3); XYZ(i,2)=(Ys1+N1(1,i)*X1(i,2)+(Ys2+N2(1,i)*X2(i,2)/2; enddisp('左片外方位元素為:Xs Ys Zs &#

15、39;);disp(WF1);disp('左片外方位元素協(xié)因素陣為：');disp(QXX1);disp('左片像點坐標(biāo)為：')disp(g1)disp('右片外方位元素為:Xs Ys Zs ');disp(WF2);disp('右片外方位元素協(xié)因素陣為：')disp(QXX2)disp('右片像點坐標(biāo)為：')disp(g2)disp('計算所得點攝影測量坐標(biāo)（X,Y,Z）為：');disp(XYZ); save 'XYZ.txt' XYZ -ascii %將計算所得結(jié)果保存到XYZ

16、.txt文件中%function Xs,Ys,Zs,q,w,k R=houfangjh(g1,Xs0,Ys0) %計算像片外方位元素%global f G m c b1 b2;f=0.152; Xs=Xs0;Ys=Ys0;Zs=m*f+G(1,3);q=0;w=0;k=0;while 1 %實現(xiàn)一個永真循環(huán)，是改正數(shù)小于限差以后跳出循環(huán) a1=cos(q)*cos(k)-sin(q)*sin(w)*sin(k); a2=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k); a3=-sin(q)*cos(w); b1_=cos(w)*sin(k); b2_=cos(w)*co

17、s(k); b3=-sin(w); c1=sin(q)*cos(k)+cos(q)*sin(w)*sin(k); c2=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k); c3=cos(q)*cos(w); R=a1,a2,a3;b1_,b2_,b3;c1,c2,c3; aX=;aY=;aZ=; for i=1:c aX(i)=a1*(G(i,1)-Xs)+b1_*(G(i,2)-Ys)+c1*(G(i,3)-Zs); aY(i)=a2*(G(i,1)-Xs)+b2_*(G(i,2)-Ys)+c2*(G(i,3)-Zs); aZ(i)=a3*(G(i,1)-Xs)+b3*

18、(G(i,2)-Ys)+c3*(G(i,3)-Zs); end xj=;yj=; for i=1:c xj(i)=-f*aX(i)/aZ(i); yj(i)=-f*aY(i)/aZ(i); end a11=;a12=;a13=;a14=;a15=;a16=; a21=;a22=;a23=;a24=;a25=;a26=; for i=1:c a11(i)=(a1*f+a3*g1(i,1)/aZ(i);a12(i)=(b1_*f+b3*g1(i,1)/aZ(i);a13(i)=(c1*f+c3*g1(i,1)/aZ(i); a21(i)=(a2*f+a3*g1(i,2)/aZ(i);a22(i)=

19、(b2_*f+b3*g1(i,2)/aZ(i);a23(i)=(c2*f+c3*g1(i,2)/aZ(i); a14(i)=g1(i,2)*sin(w)-(g1(i,1)*(g1(i,1)*cos(k)-g1(i,2)*sin(k)/f+f*cos(k)*cos(w); a15(i)=-f*sin(k)-g1(i,1)*(g1(i,1)*sin(k)+g1(i,2)*cos(k)/f; a16(i)=g1(i,2); a24(i)=-g1(i,1)*sin(w)-(g1(i,2)*(g1(i,1)*cos(k)-g1(i,2)*sin(k)/f-f*sin(k)*cos(w); a25(i)=

20、-f*cos(k)-g1(i,2)*(g1(i,1)*sin(k)+g1(i,2)*cos(k)/f; a26(i)=-g1(i,1); end lx=;ly=;for i=1:c lx(i)=g1(i,1)-xj(i);ly(i)=g1(i,2)-yj(i); end A=zeros(2*c,6); for i=1:c A(2*i-1,1)=a11(i);A(2*i-1,2)=a12(i);A(2*i-1,3)=a13(i);A(2*i-1,4)=a14(i);A(2*i-1,5)=a15(i);A(2*i-1,6)=a16(i); A(2*i,1)=a21(i); A(2*i,2)=a22

21、(i); A(2*i,3)=a23(i); A(2*i,4)=a24(i); A(2*i,5)=a25(i); A(2*i,6)=a26(i); end L=zeros(2*c,1); for i=1:c L(2*i-1,1)=lx(i); L(2*i,1)=ly(i); end X=inv(A')*A)*(A')*L; Xs=Xs+X(1,1);Ys=Ys+X(2,1);Zs=Zs+X(3,1);q=q+X(4,1);w=w+X(5,1);k=k+X(6,1); Xabs=abs(X); aaa=max(Xabs); if aaa<0.00003 %當(dāng)改正數(shù)中絕對值最大

22、的改正數(shù)小于限差0.00003 break; %后跳出循環(huán)，計算結(jié)果已經(jīng)收斂 endendglobal V;V=L'global WF QXX;WF(1)=Xs;WF(2)=Ys;WF(3)=Zs;WF(4)=q;WF(5)=w;WF(6)=k;QXX=A'*A;六 實驗結(jié)果左片外方位元素Xs,Ys,Zs, 、為：5.e+003 5.e+003 2.e+003 -7.e-005 2.e-002 9.e-002左片外方位元素協(xié)因素陣為： 4.e-008 -3.e-010 1.e-008 7.e-005 1.e-007 -2.e-006 -3.e-010 4.e-008 2.e-0

23、09 -2.e-007 7.e-005 1.e-005 1.e-008 2.e-009 1.e-008 3.e-005 6.e-006 5.e-007 7.e-005 -2.e-007 3.e-005 1.e-001 1.e-003 -1.e-003 1.e-007 7.e-005 6.e-006 1.e-003 1.e-001 3.e-002 -2.e-006 1.e-005 5.e-007 -1.e-003 3.e-002 4.e-002左片外方位元素Xs,Ys,Zs, 、為： 5.e+003 5.e+003 2.e+003 1.e-002 4.e-002 1.e-001右片外方位元素協(xié)因

24、素陣為： 3.e-008 4.e-010 -1.e-008 6.e-005 -4.e-007 1.e-006 4.e-010 3.e-008 3.e-011 -3.e-008 7.e-005 -1.e-005 -1.e-008 3.e-011 1.e-008 -2.e-005 3.e-007 -7.e-007 6.e-005 -3.e-008 -2.e-005 1.e-001 -2.e-003 4.e-003 -4.e-007 7.e-005 3.e-007 -2.e-003 1.e-001 -2.e-002 1.e-006 -1.e-005 -7.e-007 4.e-003 -2.e-002

25、 3.e-002GCP在左片和右片改正后的坐標(biāo)（x，y）為： 1.e-002 7.e-002 -7.e-002 7.e-002 8.e-002 8.e-002 -5.e-003 7.e-002 1.e-002 -7.e-002 -7.e-002 -7.e-002 8.e-002 -8.e-002 -9.e-003 -8.e-002單位權(quán)中誤差為：±1.9578e-005 所求地面點的三維坐標(biāo)（X, Y, Z）為： 5.e+003 5.e+003 5.e+002 5.e+003 5.e+003 4.e+002 5.e+003 5.e+003 5.e+002 5.e+003 5.e+00

26、3 5.e+002 5.e+003 4.e+003 4.e+002七 心得體會經(jīng)過三周的努力，這個當(dāng)初看來艱巨的任務(wù)終于在我的不懈努力下圓滿的完成了。在編程的過程中，不僅加強了我對空間后方交會求解外方位元素以及應(yīng)用空間前方交會求解物點在攝影測量坐標(biāo)系中的三維坐標(biāo)，同時也鍛煉了我們將理論知識應(yīng)用到實際中的一種能力。在這個過程中我遇到了很多的困難，經(jīng)過查看有關(guān)的資料，并在同學(xué)的幫助下成功的得到了解決。在這之前我們并沒有系統(tǒng)的學(xué)過MATLAB這個 強大的數(shù)學(xué)處理軟件，對它了解的太少，但我還是毅然的選擇了它。因為對于MATLAB，它在處理有關(guān)矩陣運算的時候顯示出了太強大的優(yōu)勢，與其他編程語言來說，簡單

27、易操作。在老師剛剛把任務(wù)交給我們的時候，我并沒與急切的就去實踐。而是將空間后方交會和空間前方交會有關(guān)的理論知識系統(tǒng)的看了幾遍。在上課的時候，老師一講好像什么都明白了似的，但是每當(dāng)合上課本的時候，好像什么都忘記了，為了能有一個比較完整的思路，我一邊看書上介紹的求解方法，一邊就把具體的求解過程用流程圖在草稿上粗略的表達了一下，這對我后面的編程過程有極大的幫助作用。在對書本知識有了一個系統(tǒng)完善的理解以后，接下了就是一個艱巨的過程了，對于理論知識我們我懂我們可以多看幾遍就可以有一個粗略的理解，但對于實踐來說，首先我們必須對理論有一個透徹的理解，然后還要將它轉(zhuǎn)化的實際應(yīng)用中去。在剛剛開始編寫代碼的時候，只是慢慢的不斷地摸索著，對于一個語句的定義，一個