单像空间后方交会实验报告(c++版).docx

单像空间后方交会姓名： 学号： 时间： Echo did this for you .2013/4/25目录一、作业任务- 3 -二、计算原理- 3 -三、算法流程- 7 -四、源程序- 8 -五、计算结果- 8 -六、结果分析- 8 -七、心得与体会- 8 -八、附页- 8 -1.c+程序- 8 -2.C+程序截图- 15 -3.matlb程序- 16 -一、 作业任务已知条件：摄影机主距f=153.24mm，x0=0，y0=0, 像片比例尺为1:40000，有四对点的像点坐标与相应的地面坐标如下表。 点号像点坐标地面坐标x(mm)y(mm)X(m)Y(m)Z(m)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31以单像空间后方交会方法，求解该像片的外方位元素。二、 计算原理1. 获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄测坐标。2. 测量控制点的像点坐标并作系统误差改正。3. 确定未知数的初始值。在竖直摄影且地面控制点大体对称分布的情况下，按如下方法确定初始值，即式中：为摄影比例尺分母；为控制点个数 4. 用三个角元素的初始值按下式计算各方向余弦值，组成旋转矩阵矩阵中各元素的计算公式如下：5. 逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标，带入以下共线方程式， 逐点计算像点坐标的近似值、6. 逐点计算误差方程式的系数和常数项，组成误差方程式。由常数项计算公式：得到常数项矩阵计算式为：7. 计算法方程的是系数矩阵和常数项，组成法方程式。8. 解法方程，求得外方位元素的改正数。9. 用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。式中：代表迭代次数。10. 将求得的外方位元素改正数与规定的限差比较，若小于限差，则迭代结束。否则用新的近似值重复49满足要求为止。11. 由误差方程的计算式：以及单位权中误差的计算式（式中：表示多余观测数）以及平差中6个参数的协因数阵最终得到参数的中误差为三、 算法流程否是否是输入原始数据像点坐标计算，系统误差改正确定外方位元素初始值组成旋转矩阵R逐点组成误差方程式并法化所有像点完否计算改正后的外方位元素计算中误差，输出成果，结束解法方程，求外方位元素改正数外方位元素改正数是否小于限差是结束并提示错误信息迭代次数小于n四、 源程序源程序代码请见附页五、 计算结果在经过三次迭代之后得到的最终成果如表1所示表1，最终计算结果Xs（米）Ys（米）Zs（米）（弧度）（弧度）（弧度）39795.45225827476.4623857572.685988-0.0039870.0021140.067578六、 结果分析由计算结果可知在拍摄照片瞬间，摄影中心在地面摄影测量坐标系中的坐标为（39795.452258，27476.462385，7572.685988）（单位：米），航向倾角为-0.003987弧度，旁向倾角为0.002114弧度，像片旋角为0.067578弧度。表2，精度评定结果参数XsYsZs中误差1.1073 1.2494 0.4881 0.0002 0.0002 0.0001 七、 心得与体会通过本次实验，我对单张像片的空间后方交会的计算原理及实现过程有了很深的认识，并在牢记各种计算公式的推导过程基础上做到了熟练应用。在本次实验中，我遇到了很多困难，但是都一一得到了解决，在解决困难的过程中的编程能力得到了提高，对其所涉及到的知识的印象也得到了加深。这个程序本身还有一些不足，例如:在迭代过程中的矩阵的相乘、输出等，没有采用编写函数的形式，而是逐一进行计算，这是程序整体看起来较臃肿。在编程过程中为了确保计算的准确性，我还利用matlab同步编写了一个相同功能的计算单张相片空间后方交会的迭代程序，并以此来验算原c+程序的正确性。现在编写完这个程序，心中很自豪，也有一些遗憾，因为我想用c#来编写，因为c#在暑期实习用过之后，时隔近一年在没有复习，不免把学会的知识又还给了老师，本想通过这次实验对c#进行复习，但是由于时间原因以及一些个人因素，没能达到自己原本的目标。总之，感谢老师给了我这次锻炼自己、提升自己的机会。八、 附页1. c+程序#includeiostreamusing namespace std;#includefstream#include iomanipconst int n=6;void inverse(double cnn);templatetypename T1,typename T2void transpose (T1*mat1,T2*mat2,int a,int b);templatetypename T1,typename T2void multi(T1*mat1,T2*mat2,T2 * result,int a,int b,int c);int main()double x42=-0.08615,-0.06899,-0.05340,0.08221,-0.01478,-0.07663,0.01046,0.06443;double X43=36589.41,25273.32,2195.17,37631.08,31324.51,728.69,39100.97,24934.98,2386.50,40426.54,30319.81,757.31;int i,j,num=0;/num为迭代次数double X06=0;/设定未知数(XS,YS,ZS,三个角度)初始值double f=0.15324;/摄影机主距f=153.24mmdouble a=1/40000.0;/像片比例尺为1:40000double R33=0;/初始化旋转矩阵R double approx_x8=0;/用于存放像点估计值。double A86=0;/定义了一个系数阵double AT68=0;/定义了A的转置矩阵double L8=0;/定义常数项const double pi=3.1415926535897932;double Asum66=0;double result26=0;double result168=0;double sumXYZ3=0;cout.precision(5);cout已知像点坐标为：n;for(i=0;i4;i+)for(j=0;j2;j+)coutfixed;if(j=0)coutxi+1= setw(8)xij ;elsecoutyi+1= setw(6)xijendl;cout已知地面四个点的坐标为：n;for(i=0;i4;i+)for(j=0;j3;j+)if(j=0)coutX;couti+1;cout= Xij ;elseif(j=1)coutY;couti+1;cout= Xij ;elsecoutZ;couti+1;cout= ;coutXijendl;coutendl;for(j=0;j3;j+)for(i=0;i4;i+)sumXYZj+=Xij;for(i=0;i2;i+)X0i=sumXYZi/4;/X0,Y0初始化X0i=1/a*f+sumXYZ2/4.0;/对Z0进行初始化doR00=cos(X03)*cos(X05)-sin(X03)*sin(X04)*sin(X05);R01=-cos(X03)*sin(X05)-sin(X03)*sin(X04)*cos(X05);R02=-sin(X03)*cos(X04);R10=cos(X04)*sin(X05);R11=cos(X04)*cos(X05);R12=-sin(X04);R20=sin(X03)*cos(X05)+cos(X03)*sin(X04)*sin(X05);R21=-sin(X03)*sin(X05)+cos(X03)*sin(X04)*cos(X05);R22=cos(X03)*cos(X04);/第一个像点的估计值，第一个点的坐标存放于X00,X01,X02approx_x0=-f*(R00*(X00-X00)+R10*(X01-X01)+R20*(X02-X02)/(R02*(X00-X00)+R12*(X01-X01)+R22*(X02-X02);approx_x1=-f*(R01*(X00-X00)+R11*(X01-X01)+R21*(X02-X02)/(R02*(X00-X00)+R12*(X01-X01)+R22*(X02-X02);/第二个像点的估计值,第2个点的坐标存放于X10,X11,X12approx_x2=-f*(R00*(X10-X00)+R10*(X11-X01)+R20*(X12-X02)/(R02*(X10-X00)+R12*(X11-X01)+R22*(X12-X02);approx_x3=-f*(R01*(X10-X00)+R11*(X11-X01)+R21*(X12-X02)/(R02*(X10-X00)+R12*(X11-X01)+R22*(X12-X02);/第三个像点的估计值,第3个点的坐标存放于X20,X21,X22approx_x4=-f*(R00*(X20-X00)+R10*(X21-X01)+R20*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02);approx_x5=-f*(R01*(X20-X00)+R11*(X21-X01)+R21*(X22-X02)/(R02*(X20-X00)+R12*(X21-X01)+R22*(X22-X02);/第四个像点的估计值,第4个点的坐标存放于X30,X31,X32approx_x6=-f*(R00*(X30-X00)+R10*(X31-X01)+R20*(X32-X02)/(R02*(X30-X00)+R12*(X31-X01)+R22*(X32-X02);approx_x7=-f*(R01*(X30-X00)+R11*(X31-X01)+R21*(X32-X02)/(R02*(X30-X00)+R12*(X31-X01)+R22*(X32-X02);for(i=0;i4;i+)/第i个像点估计值放在approx_x2*(i-1),approx_x2*i-1/*a11*/A2*i0=(R00*f+R02*approx_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a12*/A2*i1=(R10*f+R12*approx_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a13*/A2*i2=(R20*f+R22*approx_x2*i)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a21*/A2*i+10=(R01*f+R02*approx_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a22*/A2*i+11=(R11*f+R12*approx_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a23*/A2*i+12=(R21*f+R22*approx_x2*i+1)/(R02*(Xi0-X00)+R12*(Xi1-X01)+R22*(Xi2-X02);/*a14*/A2*i3=approx_x2*i+1*sin(X04)-(approx_x2*i/f*(approx_x2*i*cos(X05)-approx_x2*i+1*sin(X05)+f*cos(X05)*cos(X04);/*a15*/A2*i4=-f*sin(X05)-approx_x2*i/f*(approx_x2*i*sin(X05)+approx_x2*i+1*cos(X05);/*a16*/A2*i5=approx_x2*i+1;/*a24*/A2*i+13=-1*approx_x2*i*sin(X04)-(approx_x2*i+1/f*(approx_x2*i*cos(X05)-approx_x2*i+1*sin(X05)-f*sin(X05)*cos(X04);/*a25*/A2*i+14=-1*f*cos(X05)-approx_x2*i+1/f*(approx_x2*i*sin(X05)+approx_x2*i+1*cos(X05);/*a26*/A2*i+15=-approx_x2*i;/进行常数项的初始化for(i=0;i4;i+)L2*i=xi0-approx_x2*i;L2*i+1=xi1-approx_x2*i+1;/A的转置矩阵for(i=0;i8;i+)for(j=0;j6;j+)ATji=Aij;/实现A与AT相乘int k=0;for(i=0;i6;i+)for(j=0;j6;j+)Asumij=0;for(i=0;i6;i+)for(k=0;k6;k+)for(j=0;j8;j+)Asumik+=ATij*Ajk;/得到AT*A的逆矩阵存放在inverseAsum66中inverse(Asum);/实现矩阵Asum66与AT68的相乘,结果存放在result168中for(i=0;i6;i+)for(j=0;j8;j+)result1ij=0;for(i=0;i6;i+)for(k=0;k8;k+)for(j=0;j6;j+)result1ik+=Asumij*ATjk;/实现result168与l8的相乘，得到结果放在result26中；for(i=0;i6;i+)result2i=0;for(i=0;i6;i+)for(j=0;j8;j+)result2i+=result1ij*Lj;num+;for(i=0;i6;i+)X0i=X0i+result2i;ofstream f7(d:A.txt);f7 std:fixed;cout进行第num次迭代带得到Xs,Ys,Zs, ，改正数分别为：n;for(i=0;i6;i+)coutsetw(12)result2i;f7setw(12)result2i;coutendlendl;f7.close();getchar();while(abs(result23*206265.0)6|abs(result24*206265.0)6|abs(result25*206265.0)6);coutn满足限差条件结束循环，最终结果为：n;coutsetw(12)Xssetw(12)Yssetw(12)Zssetw(12)setw(12)setw(12)endl;ofstream f7(d:A.txt);f7 std:fixed;cout.precision(4);for(i=0;i6;i+)coutsetw(12)X0i;f7setw(16)X0i;f7.close();/今double XG61;for(i=0;i6;i+)XGi0=result2i;double AXG81,V81,VT18,VTV11,m0,D66;multi(A,XG,AXG,8,6,1); for( i=0;i8;i+) /计算改正数Vi0=AXGi0-Li; transpose (V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV00/2;coutendl;ofstream f6(d:what.txt);/f6 std:fixed;for(i=0;i6;i+)for(int j=0;j6;j+)Dij=m0*Asumij;coutsetw(10)Dij;f6setw(15)Dij;coutendl;f6endl;for(i=0;i6;i+)coutsqrt(Dii)endl;f6.close();/屏幕输出误差方程系数阵、常数项、改正数/getchar();return 0;void inverse(double cnn) int i,j,h,k;double p;double qn12;for(i=0;in;i+)/构造高斯矩阵for(j=0;jn;j+)qij=cij;for(i=0;in;i+)for(j=n;j12;j+)if(i+6=j)qij=1;elseqij=0;for(h=k=0;kn-1;k+,h+)/消去对角线以下的数据for(i=k+1;in;i+)if(qih=0)continue;p=qkh/qih;for(j=0;j12;j+) qij*=p;qij-=qkj;for(h=k=n-1;k0;k-,h-) / 消去对角线以上的数据for(i=k-1;i=0;i-)if(qih=0)continue;p=qkh/qih;for(j=0;j12;j+)qij*=p;qij-=qkj;for(i=0;in;i+)/将对角线上数据化为1p=1.0/qii;for(j=0;j12;j+)qij*=p;for(i=0;in;i+) /提取逆矩阵for(j=0;jn;j+)cij=qij+6;templatetypename T1,typename T2void transpose(T1*mat1,T2*mat2,int a,int b) int i,j; for(i=0;ib;i+) for(j=0;ja;j+)mat2ji=mat1ij;return;templatetypename T1,typename T2void multi(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c) int i,j,k;for(i=0;ia;i+)for(j=0;jc;j+)resultij=0;for(k=0;kb;k+)resultij+=mat1ik*mat2kj;return;2. C+程序截图3. matlb程序function XS,YS,ZS,phi,w,kappa=doitx_original=zeros(8,1);X=zeros(4,3);x_original=-0.08615;-0.06899;-0.05340;0.08221;-0.01478;-0.07663;0.01046;0.06443;X=36589.41,25273.32,2195.17;37631.08,31324.51,728.69;39100.97,24934.98,2386.50,;40426.54,30319.81,757.31;f=0.15324;%/摄影机主距f=153.24mma=0.000025;%/像片比例尺为1:40000m=40000;XS=38437YS=27963.155ZS=m*f+1516.9175phi=0w=0kappa=0R=zeros(3,3);num=1;dx=zeros(6,1);while num=1|abs(dx(4,1)*206265)6|abs(dx(5,1)*206265)6|abs(dx(6,1)*206265)6R(1,1)=cos(phi)*cos(kappa)-sin(phi)*sin(w)*sin(kappa);R(1,2)=-cos(phi)*sin(kappa)-sin(phi)*sin(w)*cos(kappa);R(1,3)=-sin(phi)*cos(w);R(2,1)=cos(w)*sin(kappa);R(2,2)=cos(w)*cos(kappa);R(2,3)=-sin(w);R(3,1)=sin(phi)*cos(kappa)+cos(phi)*sin(w)*sin(kappa);R(3,2)=-sin(phi)*sin(kappa)+cos(phi)*sin(w)*cos(kappa);R(3,3)=cos(phi)*cos(w);x_similar=zeros(8,1);x_similar(1,1)=-f*(R(1,1)*(X(1,1)-XS)+R(2,1)*(X(1,2)-YS)+R(3,1)*(X(1,3)-ZS)/(R(1,3)*(X(1,1)-XS)+R(2,3)*(X(1,2)-YS)+R(3,3)*(X(1,3)-ZS);%x1x_similar(2,1)=-f*(R(1,2)*(X(1,1)-XS)+R(2,2)*(X(1,2)-YS)+R(3,2)*(X(1,3)-ZS)/(R(1,3)*(X(1,1)-XS)+R(2,3)*(X(1,2)-YS)+R(3,3)*(X(1,3)-ZS);%y1x_similar(3,1)=-f*(R(1,1)*(X(2,1)-XS)+R(2,1)*(X(2,2)-YS)+R(3,1)*(X(2,3)-ZS)/(R(1,3)*(X(2,1)-XS)+R(2,3)*(X(2,2)-YS)+R(3,3)*(X(2,3)-ZS);%x2x_similar(4,1)=-f*(R(1,2)*(X(2,1)-XS)+R(2,2)*(X(2,2)-YS)+R(3,2)*(X(2,3)-ZS)/(R(1,3)*(X(2,1)-XS)+R(2,3)*(X(2,2)-YS)+R(3,3)*(X(2,3)-ZS);%y2x_similar(5,1)=-f*(R(1,1)*(X(3,1)-XS)+R(2,1)*(X(3,2)-YS)+R(3,1)*(X(3,3)-ZS)/(R(1,3)*(X(3,1)-XS)+R(2,3)*(X(3,2)-YS)+R(3,3)*(X(3,3)-ZS);%x3x_similar(6,1)=-f*(R(1,2)*(X(3,1)-XS)+R(2,2)*(X(3,2)-YS)+R(3,2)*(X(3,3)-ZS)/(R(1,3)*(X(3,1)-XS)+R(2,3)*(X(3,2)-YS)+R(3,3)*(X(3,3)-ZS);%y3x_similar(7,1)=-f*(R(1,1)*(X(4,1)-XS)+R(2,1)*(X(4,2)-YS)+R(3,1)*(X(4,3)-ZS)/(R(1,3)*(X(4,1)-XS)+R(2,3)*(X(4,2)-YS)+R(3,3)*(X(4,3)-ZS);%x4x_similar(8,1)=-f*(R(1,2)*(X(4,1)-XS)+R(2,2)*(X(4,2)-YS)+R(3,2)*(X(4,3)-ZS)/(R(1,3)*(X(4,1)-XS)+R(2,3)*(X(4,2)-YS)+R(3,3)*(X(4,3)-ZS);%y4a=zeros(8,6);for i=1:4 Z=(R(1,3)*(X(i,1)-XS)+R(2,3)*(X(i,2)-YS)+R(3,3)*(X(i,3)-ZS); a(i*2-1,1)=1/Z*(R(1,1)*f+R(1,3)*x_sim