摄影测量-空间前交、后交

1、- 前交程序设计（实验报告 ）姓名：班级：学号：时间：空间后交-前交程序设计一、实验目的用C 、VB或MATLAB语言编写空间后方交会-空间前方交会程序提交实习报告：程序框图、程序源代码、计算结果、体会计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度、实验数据点号左片右片地面摄影测量坐标XyXyX1GCP116.01279. 963-7X9378.7065083. 2055852.099527.925GCP288. 5681.154-5.25278. 1845780. 025906.36557L 549GCP31X362-79.37-79. 1221-78.8795210.879

2、4258. 44646L81GCP482. 2480. 027-9.8871-80. 0895909. 2644314, 283455. 484151.75880. 555-39.95378. 463214 618-0.231-76. 0060. Q36349. 88-0.7EJ2-42. 2。1-1. 022486. 14L 34£-7. 706-2.112548. 035-79.962-44. 438-79.736f=150.000mm, x0=0, y0=0三、实验思路1利用空间后方交会求左右像片的外方位元素获取m(于像片中选取两点，于地面摄影测量坐标系中选取同点,分别计算距离

3、，距离比值即为 m), x,y,f,X,Y,Z(2)确定未知数初始值 Xs,Ys,Zs,q,w,k计算旋车矩阵R(4)逐点计算像点坐标的近似值(x) , (y)(5)组成误差方程式(6)组成法方程式(7)解求外方位元素(8)检查是否收敛，即将求得的外方位元素的改正数与规定限差比较, 小于限差即终止；否则用新的近似值重复步骤(3)-(7)2利用求出的外方位元素进行空间前交，求出待定点地面坐标(1)用各自像片的角元素计算出左、右像片的方向余弦值，组成旋转矩阵R1,R2Bz(2)根据左、右像片的外方位元素，计算摄影基线分量Bx, By,计算像点的像空间辅助坐标(X1,Y1,Z1)和(X2,Y2,Z2

4、)(4)计算点投影系数N1和N2(5)计算未知点的地面摄影测量坐标四、实验过程程序框图testvarAAndL输入/ spacehoujiao* /计算spaceqianiao输出deg2dmsxyok程序代码函数AandL精选范本%求间接平差时需要的系数% 已知%a=像点坐标*扫=像点坐标y, f内方位元素主距%j =q, ° =w, k =k%像空间坐标系X,Y,Z%地面摄影测量坐标系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

5、)*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);% 共线方程的分子分母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

6、-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+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(

7、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;% 组成一个矩阵，并返回A1=a11,a12,a13,a14,a15,a16;A2=a21,a22,a23,a24,a25,a26;L1=lx;L2=ly;函数 deg2dms% 角度转度分秒function y=deg2dms(x) a=floor(x);b=floor(x-a)*60);c=(x-a-b/60)*3600;y=a+(b/100)+(c/10000);函数 dms2deg%度分秒转度function y=d

8、ms2deg(x) a=floor(x);b=floor(x-a)*100);c=(x-a-b/100)*10000;y=a+b/60+c/3600;函数 ok% 目的是为了保证各取的值的有效值%xy 为 n*1,a 为 1*n function result=ok(xy,a) format short gi=size(xy,1);for n=1:io=xy(n)-floor(xy(n,1);o=round(o*(10Aa(n)/(10Aa(n);xy(n,1)=floor(xy(n,1)+o;endformat long g result=xy;函数 rad2dmsxy%求度分秒表现形式的三

9、个外方位元素，三个角度function xydms=rad2dmsxy(xy) a,b,c,d,e,f=testvar(xy);d=deg2dms(rad2deg(d);e=deg2dms(rad2deg(e);f=deg2dms(rad2deg(f); xydms=a,b,c,d,e,f'函数 spacehoujiao% 空间后交% f精选范本%输入 p( 2*n,1)%像点坐标x,y,X,Y,Z,均为(n,1)function xy,m,R=spacehoujiao(p,x,y,f,X,Y,Z) format long;%权的矢量化 ,这是等精度时的，如果非，将函数参数改为 PP=

10、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)42+(y(2)-y(1)A2);% 两像点之间距离Sd=sqrt(X(2)-X(1)A2+(Y(2)-Y42);%两地面控制点之间距离m=Sd/Sx; % 图像比例系数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),

11、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;% 循环体while 1%dXs,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&l

12、t;0.01)&&(detq<pi/648000)&&(detw<pi/648000)&&(detk<pi/648000)break;elseV=(a*det-l);Q=inv(a'*P*a);m0=m0+sqrt(V'*P*V)/(2*j-6);%m0 需要每次的改正数算出来相加 %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)

13、,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);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

14、=q+dq;w=w+dw;k=k+dk;% 可以输出迭代次数的 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),mo(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);%旋转矩阵函数 space

15、qianjiao%空间前交%输入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(x1,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),

16、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(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函数 testvar%分割矩阵。%将矩阵的每行元素打包给元素。%用法 Xs1,Ys1,Zs1,q1,w1,k1=testvar(xy1);function vara

17、rgout=testvar(arrayin)for k=1:nargoutvarargoutk=arrayin(k,:);end函数 xyR%计算旋转矩阵，通过六个内方位元素%xy ( 6,1)function R=xyR(xy)a,b,c,q,w,k=testvar(xy);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)+

18、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;%空间后交%像片坐标x11=(1e-3)*16.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;%地面摄影测量坐标X1=5083

19、.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六个内方位元素矩阵，xy1,左片xy2,右片%m1,m2误差矩阵，ml,左片m2,右片%R1,R2旋转矢I阵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(rad2dm

20、sxy(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,R

21、2,xy1,xy2);X2d=(ok(X2',3,3,3,3,3)'Y2d=(ok(Y2',3,3,3,3,3)'Z2d=(ok(Z2',3,3,3,3,3)'(3)计算结果左片六个外方位元素xyldmsJI位中误差和六个外方位元素中误差 ml左片六个外方位元素xy2dm浦位中误差和六个外方位元素中误差 m2空间前交得到的坐标X2d, Y2dz2dsyldns =4999. U4999.7292000. 0020. 00441. 39555. ”26Ey2djus =5896.8295 DTD. 2442030.430.49372.3312仅 1

22、946Q. 00046317741663986 IT. 23522S130S27S 22. 44.00622210581 7.55104111564557 0.0109507237615990 0.010B12991M90D95 0.00409513286115f3哈=口，OlDOfi6M2296768690125, 681444013600L33. 61 29779"" 9_475890532378?20. 0100475421635183 0.016397-30 526718510.CD548053994E78844堂C -6491.8051d7. SBBT2rl -一

23、讯 口08国眨55时d =MX 7234组 DDL4D5. 787555LM1082.7275109. B02"轨 ias677523. 424(3,623点号左片右声地面摄影测量坐标.XyyXYZGCP116. 012Td 963-73. 93TS. 7OG5OS3- 2055852_ 8952T. 90SHK 56SI- ±34-h.V«L ±«40239g. MSbh71- b4<3GCT313.362Td 37-79l 122-7«, 879521UL 6794258. 14546K 811 GCF424SBD. rarF,