攝影測(cè)量-空間前交、后交

1、空間后交-前交程序設(shè)計(jì)(實(shí)驗(yàn)報(bào)告)姓名：班級(jí)：學(xué)號(hào)：時(shí)間：空間后交-前交程序設(shè)計(jì)一、實(shí)驗(yàn)?zāi)康挠?C 、VB或MATLAB語(yǔ)言編寫(xiě)空間后方交會(huì)-空間前方交會(huì)程序提交實(shí)習(xí)報(bào)告：程序框圖、程序源代碼、計(jì)算結(jié)果、體會(huì)計(jì)算結(jié)果：像點(diǎn)坐標(biāo)、地面坐標(biāo)、單位權(quán)中誤差、外方位元素及其精度二、實(shí)驗(yàn)數(shù)據(jù)f=150.000mm，x0=0，y0=0三、實(shí)驗(yàn)思路1.利用空間后方交會(huì)求左右像片的外方位元素(1).獲取m(于像片中選取兩點(diǎn)，于地面攝影測(cè)量坐標(biāo)系中選取同點(diǎn)，分別計(jì)算距離，距離比值即為m)，x,y,f,X,Y,Z(2).確定未知數(shù)初始值Xs,Ys,Zs,q,w,k(3).計(jì)算旋轉(zhuǎn)矩陣R(4).逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近

2、似值（x），（y）(5).組成誤差方程式(6).組成法方程式(7).解求外方位元素(8).檢查是否收斂，即將求得的外方位元素的改正數(shù)與規(guī)定限差比較，小于限差即終止；否則用新的近似值重復(fù)步驟（3）-（7）2.利用求出的外方位元素進(jìn)行空間前交，求出待定點(diǎn)地面坐標(biāo)(1).用各自像片的角元素計(jì)算出左、右像片的方向余弦值，組成旋轉(zhuǎn)矩陣R1,R2(2).根據(jù)左、右像片的外方位元素 ，計(jì)算攝影基線分量Bx，By，Bz(3).計(jì)算像點(diǎn)的像空間輔助坐標(biāo)（X1,Y1,Z1）和（X2,Y2,Z2）(4).計(jì)算點(diǎn)投影系數(shù)N1和N2(5).計(jì)算未知點(diǎn)的地面攝影測(cè)量坐標(biāo)四、實(shí)驗(yàn)過(guò)程程序框圖testvar輸入AAndLsp

3、acehoujiao計(jì)算spaceqianiaodeg2dmsxy輸出ok程序代碼函數(shù)AandL%求間接平差時(shí)需要的系數(shù)%已知%a=像點(diǎn)坐標(biāo)x,b=像點(diǎn)坐標(biāo)y，f內(nèi)方位元素主距%=q,=w,=k%像空間坐標(biāo)系X,Y,Z%地面攝影測(cè)量坐標(biāo)系Xs,Ys,Zsfunction A1,L1,A2,L2=AandL(a,b,f,q,w,k,X,Y,Z,Xs,Ys,Zs)%選擇矩陣元素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);

4、b2=cos(w)*cos(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);%共線方程的分子分母X_=a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs);Y_=a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs);Z_=a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs);%近似值x=-f*X_/Z_;y=-f*Y_/Z_;%A組成L組成a11=1/Z_*(a1*f+a3*x);a12=1/Z_*(b1*f+

5、b3*x);a13=1/Z_*(c1*f+c3*x);a21=1/Z_*(a2*f+a3*y);a22=1/Z_*(b2*f+b3*y);a23=1/Z_*(c2*f+c3*y);a14=y*sin(w)-(x/f*(x*cos(k)-y*sin(k)+f*cos(k)*cos(w);a15=-f*sin(k)-x/f*(x*sin(k)+y*cos(k); a16=y;a24=-x*sin(w)-(y/f*(x*cos(k)-y*sin(k)-f*sin(k)*cos(w);a25=-f*cos(k)-y/f*(x*sin(k)+y*cos(k); a26=-x;lx=a-x;ly=b-y;

6、%組成一個(gè)矩陣，并返回A1=a11,a12,a13,a14,a15,a16;A2=a21,a22,a23,a24,a25,a26;L1=lx;L2=ly;函數(shù)deg2dms%角度轉(zhuǎn)度分秒function y=deg2dms(x)a=floor(x);b=floor(x-a)*60);c=(x-a-b/60)*3600;y=a+(b/100)+(c/10000);函數(shù)dms2deg%度分秒轉(zhuǎn)度f(wàn)unction y=dms2deg(x)a=floor(x);b=floor(x-a)*100);c=(x-a-b/100)*10000;y=a+b/60+c/3600;函數(shù)ok%目的是為了保證各取的值的

7、有效值%xy為n*1,a為1*nfunction result=ok(xy,a)format short gi=size(xy,1); for n=1:io=xy(n)-floor(xy(n,1);o=round(o*(10a(n)/(10a(n);xy(n,1)=floor(xy(n,1)+o;endformat long gresult=xy;函數(shù)rad2dmsxy%求度分秒表現(xiàn)形式的三個(gè)外方位元素，三個(gè)角度f(wàn)unction xydms=rad2dmsxy(xy)a,b,c,d,e,f=testvar(xy);d=deg2dms(rad2deg(d);e=deg2dms(rad2deg(e

8、);f=deg2dms(rad2deg(f);xydms=a,b,c,d,e,f'函數(shù)spacehoujiao%空間后交% f%輸入p（2*n,1）%像點(diǎn)坐標(biāo)x,y,X,Y,Z,均為（n,1）function xy,m,R=spacehoujiao(p,x,y,f,X,Y,Z)format long; %權(quán)的矢量化,這是等精度時(shí)的，如果非，將函數(shù)參數(shù)改為PP=diag(p);%求nj=size(X,2);%初始化Xs=0;Ys=0;Zs=0;for n=1:jXs=Xs+X(n);Ys=Ys+Y(n);Zs=Zs+Z(n);endSx=sqrt(x(2)-x(1)2+(y(2)-y(1

9、)2);%兩像點(diǎn)之間距離Sd=sqrt(X(2)-X(1)2+(Y(2)-Y(1)2);%兩地面控制點(diǎn)之間距離m=Sd/Sx; %圖像比例系數(shù)Xs=Xs/j;Ys=Ys/j;Zs=m*f+Zs/j;m0=0;q=0;w=0;k=0;i=0;a=rand(2*j,6);l=rand(2*j,1);%for n=1:ja(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n),Xs,Ys,Zs);enddet=inv(a'*P*a)*transpose(a)*P*l;%循環(huán)體while 1%d

10、Xs,dYs,dZs,dq,dw,dk=testvar(det);detXs=abs(dXs);detYs=abs(dYs);detZs=abs(dZs);detq=abs(dq);detw=abs(dw);detk=abs(dk);%if (detXs<0.01)&&(detYs<0.01)&&(detZs<0.01)&&(detq<pi/648000)&&(detw<pi/648000)&&(detk<pi/648000)break;else V=(a*det-l);Q=in

11、v(a'*P*a);m0=m0+sqrt(V'*P*V)/(2*j-6);%m0需要每次的改正數(shù)算出來(lái)相加%Xs=Xs+dXs;Ys=Ys+dYs;Zs=Zs+dZs;q=q+dq;w=w+dw;k=k+dk;%for n=1:ja(2*n-1,:),l(2*n-1,1),a(2*n,:),l(2*n,1)=AandL(x(n),y(n),f,q,w,k,X(n),Y(n),Z(n),Xs,Ys,Zs);enddet=inv(a'*P*a)*transpose(a)*P*l;i=i+1;%end%enddXs,dYs,dZs,dq,dw,dk=testvar(det);

12、detXs=abs(dXs);detYs=abs(dYs);detZs=abs(dZs);detq=abs(dq);detw=abs(dw);detk=abs(dk);V=(a*det-l);Q=inv(a'*P*a);m0=m0+sqrt(V'*P*V)/(2*n-6);%Xs=Xs+dXs;Ys=Ys+dYs;Zs=Zs+dZs;q=q+dq;w=w+dw;k=k+dk;%可以輸出迭代次數(shù)的i%Xs,Ys,Zs,q,w,k,i,dXs,dYs,dZs,dq,dw,dk,detXs,detYs,detZs%精度mo=m0*sqrt(Q);m=mo(1,1),mo(2,2),m

13、o(3,3),mo(4,4),mo(5,5),mo(6,6)'mXs,mYs,mZs,mq,mw,mk=testvar(m);%輸出xy=Xs,Ys,Zs,q,w,k'%輸出（6,1）的外方位元素m=m0,mXs,mYs,mZs,mq,mw,mk'%單位誤差，各元素中誤差R=xyR(xy);%旋轉(zhuǎn)矩陣函數(shù)spaceqianjiao%空間前交%輸入f%輸入x1,y1,x2,y2,R1,R2,xy1,xy2 （n,1）%輸出X,Y,Z （n,1）function X,Y,Z=spaceqianjiao(x1,y1,x2,y2,f,R1,R2,xy1,xy2)i=size(x

14、1,2);Xs1,Ys1,Zs1,q1,w1,k1=testvar(xy1);Xs2,Ys2,Zs2,q2,w2,k2=testvar(xy2);for n=1:iX1(n),Y1(n),Z1(n)=testvar(R1*x1(n),y1(n),-f');X2(n),Y2(n),Z2(n)=testvar(R2*x2(n),y2(n),-f');Bx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1;N1=(Bx*Z2(n)-Bz*X2(n)/(X1(n)*Z2(n)-X2(n)*Z1(n);N2=(Bx*Z1(n)-Bz*X1(n)/(X1(n)*Z2(n)-X2(

15、n)*Z1(n);X(n)=Xs1+N1*X1(n);Z(n)=Zs1+N1*Z1(n);Y(n)=0.5*(Ys1+N1*Y1(n)+(Ys2+N2*Y2(n);end函數(shù)testvar%分割矩陣。%將矩陣的每行元素打包給元素。%用法Xs1,Ys1,Zs1,q1,w1,k1=testvar(xy1);function varargout=testvar(arrayin)for k=1:nargout varargoutk=arrayin(k,:);end函數(shù)xyR%計(jì)算旋轉(zhuǎn)矩陣，通過(guò)六個(gè)內(nèi)方位元素%xy （6,1）function R=xyR(xy)a,b,c,q,w,k=testvar(x

16、y);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)*cos(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;命令行f=0.15;%空間后交%像片坐標(biāo)x11=(1e-3)*16.

17、012,88.56,13.362,82.24;y11=(1e-3)*79.963,81.134,-79.37,-80.027;x21=(1e-3)*-73.93,-5.252,-79.122,-9.887;y21=(1e-3)*78.706,78.184,-78.879,-80.089;%地面攝影測(cè)量坐標(biāo)X1=5083.205,5780.02,5210.879,5909.264;Y1=5852.099,5906.365,4258.446,4314.283;Z1=527.925,571.549,461.81,455.484;p=1,1,1,1,1,1,1,1;%xy1,xy2 六個(gè)內(nèi)方位元素矩陣

18、，xy1,左片xy2，右片%m1,m2 誤差矩陣，m1,左片m2，右片%R1,R2 旋轉(zhuǎn)矩陣R1,左片R2，右片xy1,m1,R1=spacehoujiao(p,x11,y11,f,X1,Y1,Z1);xy2,m2,R2=spacehoujiao(p,x21,y21,f,X1,Y1,Z1);xy1dms=ok(rad2dmsxy(xy1),3,3,3,4,4,4)m1xy2dms=ok(rad2dmsxy(xy2),3,3,3,4,4,4)m2%空間前交x12=(1e-3)*51.758,14.618,49.88,86.14,48.035;y12=(1e-3)*80.555,-0.231,-0.782,-1.346,-79.962;x22=(1e-3)*-39.953,-76.006,-42.201,-7.706,-44.438;y22=(1e-3)*78.463,0.036,-1.022,-2.112,-79.736;X2,Y2,Z2=spaceqianjiao(x12,y12,x22,y22,f,R1,R2,xy1,xy2);X2d=(ok(X2',3,3,3,3,3)'Y2d=(ok(Y2',