《摄 影 测 量 学》单像空间后方交会实习报告 班 级： XXXX00000 姓 名： X X X000 学 号：XXXXXXXXXXXXX 指导教师： X X X 一、 实习目的 1、掌握空间后方交会的定义和实现算法； 2、了解摄影测量平差的基本过程； 3、熟练MATLAB等程序编写二、 实习原理利用至少三个已知地面控制点的坐标、、，与其影像上对应的三个像点的影像坐标、、，根据共线方程，反求该像片的外方位元素共线条件方程式： 将共线式线性化并取一次小值项得：三、 解算过程①获取已知数据包括影像比例尺1/m，平均摄影距离（航空摄影的航高）H，内方位元素x0、y0、f，控制点的空间坐标X、Y、Z ②量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标③确定未知数的初始值单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，Xs0和Ys0为均值，Zs0为航高，φ、ω、κ的初值都设为0。
或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出 ④计算旋转矩阵R利用角元素近似值计算方向余弦值，组成R阵⑤逐点计算像点坐标的近似值利用未知数的近似值按共线条件式计算控制点像点坐标的近似值（x），(y)⑥逐点计算误差方程式的系数和常数项，组成误差方程式⑦计算法方程的系数矩阵ATA与常数项ATL，组成法方程式⑧解求外方位元素根据法方程，解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值 ⑨检查计算是否收敛将所求得的外方位元素的改正数与规定的限差比较，通常对φ，ω，κ的改正数△φ，△ω，△κ给予限差，通常为0.000001弧度，当3个改正数均小于0.000001弧度时，迭代结束否则用新的近似值重复④～⑧步骤的计算，直到满足要求为止四、程序框图五、主要代码[fn,pn,fi]=uigetfile('*.txt','选择文件');I=load([pn fn]);x=I(:,1)';y=I(:,2)';X=I(:,3)';Y=I(:,4)';Z=I(:,5)';s=0;S=0;for i=1:3; j=i+1; sij=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2); Sij=sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2); s=s+sij; S=S+Sij;endm=S*1000/s;f=153.24;Xs0=(X(1)+X(2)+X(3)+X(4))/4;Ys0=(Y(1)+Y(2)+Y(3)+Y(4))/4;Zs0=m*f;fai0=0;omig0=0;ka0=0;for v=1:+inf; a1=cos(fai0)*cos(ka0)-sin(fai0)*sin(omig0)*sin(ka0); a2=-cos(fai0)*sin(ka0)-sin(fai0)*sin(omig0)*cos(ka0); a3=-sin(fai0)*cos(omig0); b1=cos(omig0)*sin(ka0); b2=cos(omig0)*cos(ka0); b3=-sin(omig0); c1=sin(fai0)*cos(ka0)+cos(fai0)*sin(omig0)*sin(ka0); c2=-sin(fai0)*sin(ka0)+cos(fai0)*sin(omig0)*cos(ka0); c3=cos(fai0)*cos(omig0); R=[a1 a2 a3;b1 b2 b3;c1 c2 c3]; l=[];A=[]; for h=1:4; O=a1*(X(h)-Xs0)+b1*(Y(h)-Ys0)+c1*(Z(h)-Zs0); P=a2*(X(h)-Xs0)+b2*(Y(h)-Ys0)+c2*(Z(h)-Zs0); Q=a3*(X(h)-Xs0)+b3*(Y(h)-Ys0)+c3*(Z(h)-Zs0); x1=-f*O/Q; y1=-f*P/Q; a11=(a1*f+a3*x(h))/Q; a12=(b1*f+b3*x(h))/Q; a13=(c1*f+c3*x(h))/Q; a14=y(h)*sin(omig0)-(x(h)/f*(x(h)*cos(ka0)-y(h)*sin(ka0))+f*cos(ka0))*cos(omig0); a15=-f*sin(ka0)-x(h)/f*(x(h)*sin(ka0)+y(h)*cos(ka0)); a16=y(h); a21=(a2*f+a3*y(h))/Q; a22=(b2*f+b3*y(h))/Q; a23=(c2*f+c3*y(h))/Q; a24=-x(h)*sin(omig0)-(y(h)/f*(x(h)*cos(ka0)-y(h)*sin(ka0))-f*sin(ka0))*cos(omig0); a25=-f*cos(ka0)-y(h)/f*(x(h)*sin(ka0)+y(h)*cos(ka0)); a26=-x(h); lx=x(h)-x1; ly=y(h)-y1; lh=[lx ly]'; Ah=[a11,a12,a13,a14,a15,a16;a21,a22,a23,a24,a25,a26]; A=[A;Ah]; l=[l;lh]; end XX=inv(A'*A)*A'*l; Xs0=Xs0+XX(1); Ys0=Ys0+XX(2); Zs0=Zs0+XX(3); fai0=fai0+XX(4); omig0=omig0+XX(5); ka0=ka0+XX(6); R=[a1,a2,a3;b1,b2,b3;c1,c2,c3] Xs=Xs0 Ys=Ys0 Zs=Zs0 fai=fai0 omig=omig0 ka=ka0 If abs(XX(4))<0.0000291&&abs(XX(5))<0.0000291&&abs(XX(6))<0.0000291 break endEndR=[a1,a2,a3;b1,b2,b3;c1,c2,c3];V=A*XX-l;Qii=inv(A'*A);m0=sqrt(V'*V/2);mi=m0*sqrt(Qii);m=diag(mi);m=[m(1),m(2),m(3),m(4)/pi*180*3600,m(5)/pi*180*3600,m(6)/pi*180*3600];mXs=m(1);mYs=m(2);mZs=m(3);mfai=m(4);momig=m(5);mka=m(6);fp=fopen('后方交会计算结果.txt','wt'); fprintf(fp,'迭代次数:%d\n',v); fprintf(fp,'\n旋转矩阵R：\n');[m,n]=size(R); for i=1:1:m for j=1:1:n if j==n fprintf(fp,'%g\n',R(i,j)); else fprintf(fp,'%g\t',R(i,j)); end endendfprintf(fp,'\n外方位元素解：\n');fprintf(fp,'Xs=%g\t',Xs);fprintf(fp,'Ys=%g\t',Ys);fprintf(fp,'Zs=%g\n',Zs);fprintf(fp,'φ=%g\t',fai);fprintf(fp,'ω=%g\t',omig);fprintf(fp,'κ=%g\n',ka);fprintf(fp,'\n单位权中误差的值：mo=%g\n',m0);fprintf(fp,'\n外方位元素中误差为：\n');fprintf(fp,'mXs=%g米\t',mXs);fprintf(fp,'mYs=%g米\t',mYs);fprintf(fp,'mZs=%g米\n',mZs);fprintf(fp,'mφ=%g秒\t',mfai);fprintf(fp,'mω=%g秒\t',momig);fprintf(fp,'mκ=%g秒\n',mka);六、 实习结果 1、起算数据 2、运行结果 3、几次迭代的数据。