摄影测量学单像空间后方交会编程实习报告

摄影测量学单像空间后方交会编程实习报告班级：130x姓名：xx学号：2013302590xxx指导老师：李欣一、实习目的通过对提供的数据进行计算，输出像片的外方位元素并评定精度。深入理解单像空间后方交会的思想，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过尝试编程实现加强编程处理问题的能力和对实习内容的理解，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。了解摄影测量平差的基本过程，掌握空间后方交会的定义和实现算法。二、实习内容根据学习的单像空间后方交会的知识，用程序设计语言（C++或C语言）编写一个完整的单像空间后方交会程序，通过对提供的数据进行计算，输出像片的外方位元素并评定精度。三、实习数据已知航摄仪的内方位元素：fk=153.24mm，x0=y0=0，摄影比例尺为1:15000；4个地面控制点的地面坐标及其对应像点的像片坐标：像片坐标地面点坐标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四、实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。五、实习流程1.获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，𝑓；获取控制点的空间坐标Xt，Yt，Zt。2.量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。3.确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：𝑍𝑠0=H=m∙𝑓𝑋𝑠0=1𝑛∑𝑋𝑡𝑖𝑛𝑖−1𝑌𝑠0=1𝑛∑𝑌𝑡𝑖𝑛𝑖−1𝜑0=𝜔0=0式中：m为摄影比例尺分母,n为控制点个数；4.计算旋转矩阵R。利用角元素的近似值按公式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。5.逐点计算像点坐标的近似值。利用未知数的近似值按共线条件方程计算控制点像点坐标的近似值(x),(y)。6.逐点计算误差方程式的系数和常数项，组成误差方程。7.计算法方程的系数矩阵ATA与常数项ATL，组成法方程。8.解求外方位元素。根据法方程，解求外方位元素的改正数，并与相应的近似值求和，得到外方位元素新的近似值。9.检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对ϕ,ω,κ的改正数△φ,△ω,△κ给予限差，当改正数小于限差时，迭代结束。否则用新的近似值重复（4）——（8）步骤计算，直到满足要求为止。10.空间后方交会的精度估计：按上述方法所求得的影像外方位元素的精度可以通过法方程式中未知数的系数矩阵的逆阵（ATA）-1来解求，此时视像点坐标为等精度不相关观测值。因为(ATA）-1中第i个主对角线上的元素Qii就是法方程式中第i个未知数的权倒数，若单位权中误差为m0，则第i个未知数的中误差为：𝑚𝑖=√𝑄𝑖𝑖𝑚0当参加空间后方交会的控制点有n个时，则单位权中误差可按下式计算：𝑚0=±√[𝑣𝑣]2𝑛−6流程图如下:六、程序代码#includestdio.h#includeMatrix.h//矩阵运算头文件，来自网络，已上传voidmain(){doubleXs,Ys,Zs,q,w,k;//外方位元素doublex0,y0,f;//内方位元素doublex[4],y[4];//影像坐标doubleX0[4],Y0[4],Z0[4];//地面坐标doublem;//比例尺doublea1,a2,a3,b1,b2,b3,c1,c2,c3;//旋转矩阵intn=0;//迭代次数//地面坐标X0[0]=36589.41;X0[1]=37631.08;X0[2]=39100.97;X0[3]=40426.54;Y0[0]=25273.32;Y0[1]=31324.51;Y0[2]=24934.98;Y0[3]=30319.81;Z0[0]=2195.17;Z0[1]=728.69;Z0[2]=2386.50;Z0[3]=757.31;//影像坐标x[0]=-0.08615;x[1]=-0.05340;x[2]=-0.01478;x[3]=0.01046;y[0]=-0.06899;y[1]=0.08221;y[2]=-0.07663;y[3]=0.06443;//确定内外方位元素初始值x0=0;y0=0;f=153.24/1000;m=15000;Xs=0;Ys=0;q=0;w=0;k=0;for(inti=0;i4;i++){Xs+=X0[i];Ys+=Y0[i];}Xs/=4;Ys/=4;Zs=f*m;MatrixA(8,6);MatrixL(8,1);MatrixAT(6,8);MatrixATA(6,6);MatrixATA_(6,6);MatrixXx(6,1);//迭代计算do{//旋转矩阵Ra1=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);//计算像点坐标for(inti=0;i4;i++){doubleX=a1*(X0[i]-Xs)+b1*(Y0[i]-Ys)+c1*(Z0[i]-Zs);doubleY=a2*(X0[i]-Xs)+b2*(Y0[i]-Ys)+c2*(Z0[i]-Zs);doubleZ=a3*(X0[i]-Xs)+b3*(Y0[i]-Ys)+c3*(Z0[i]-Zs);doublexx=x[i]-x0;doubleyy=y[i]-y0;A[2*i][0]=(a1*f+a3*xx)/Z;A[2*i][1]=(b1*f+b3*xx)/Z;A[2*i][2]=(c1*f+c3*xx)/Z;A[2*i][3]=yy*sin(w)-(xx*(xx*cos(k)-yy*sin(k))/f+f*cos(k))*cos(w);A[2*i][4]=-f*sin(k)-xx*(xx*sin(k)+yy*cos(k))/f;A[2*i][5]=yy;A[2*i+1][0]=(a2*f+a3*yy)/Z;A[2*i+1][1]=(b2*f+b3*yy)/Z;A[2*i+1][2]=(c2*f+c3*yy)/Z;A[2*i+1][3]=-xx*sin(w)-(yy*(xx*cos(k)-yy*sin(k))/f-f*sin(k))*cos(w);A[2*i+1][4]=-f*cos(k)-yy*(xx*sin(k)+yy*cos(k))/f;A[2*i+1][5]=-xx;L[2*i][0]=x[i]-(-f*X/Z);L[2*i+1][0]=y[i]-(-f*Y/Z);}//法方程的解X=(ATA)_ATLAT=A.Transpose();//转置ATA=AT*A;ATA_=ATA.Inverse();//求逆Xx=(ATA_*AT)*L;//法方程解Xs+=Xx[0][0];Ys+=Xx[1][0];Zs+=Xx[2][0];q+=Xx[3][0];w+=Xx[4][0];k+=Xx[5][0];n++;//printf(%d\n,n);//printf(%lf%lf%lf%.5lf%.5lf%.5lf\n,Xs,Ys,Zs,q,w,k);}while((abs(Xx[3][0])0.000001||abs(Xx[4][0])0.000001||abs(Xx[5][0])0.000001)&&n100);//求中误差doublem0,mi[6];//单位权中误差的绝对值m0以及第i个未知数的中误差mim0=sqrt(((A*Xx-L).Transpose()*(A*Xx-L)).ToDouble()/(2*4-6));MatrixQ(6,6);Q=(A.Transpose()*A).Inverse();for(inti=0;i6;i++){mi[i]=(sqrt(Q[i][i]))*m0;}//printf(%lf\n,m0);//输出printf(迭代次数:%d\n,n);printf(旋转矩阵R:\n);printf(%.5lf%.5lf%.5lf\n,a1,a2,a3);printf(%.5lf%.5lf%.5lf\n,b1,b2,b3);printf(%.5lf%.5lf%.5lf\n,c1,c2,c3);printf(\n外方位元素解:\nWs=%.2lfYs=%.2lfZs=%.2lf\nq=%.5lfw=%.5lfk=%.5lf\n\n,Xs,Ys,Zs,q,w,k);printf(单位权中误差的绝对值:%lfm\n,m0);printf(Xs的精度:%lfm\n,mi[0]);printf(Ys的精度:%lfm\n,mi[1]);printf(Zs的精度:%lfm\n,mi[2]);printf(q的精度:%lf\n,mi[3]);printf(w的精度:%lf\n,mi[4]);printf(k的精度:%lf\n,mi[5]);//保存到文件，结果文件默认保存在程序目录FILE*fp;fp=fopen(result.txt,w);//fprintf(fp,迭代次数:%d\n,n);fprintf(fp,旋转矩阵R:\n);fprintf(fp,%.5lf%.5lf%.5lf\n,a1,a2,a3);fprintf(fp,%.5lf%.5lf%.5lf\n,b1,b2,b3);fprintf(fp,%.5lf%.5lf%.5lf\n,c1,c2,c3);fprintf(fp,\n外方位元素解:\nWs=%.2lfYs=%.2lfZs=%.2lf\nq=%.5lfw=%.5lfk=%.5lf\n\n,Xs,Ys,Zs,q,w,k);fprintf(fp,单位权中误差的绝对值:%lfm\n,m0);fprintf(fp,Xs的精度:%lfm\n,mi[0]);fprintf(fp,Ys的精度:%lfm\n,mi[1]);fprintf(fp,Zs的精度:%lfm\n,mi[2]);fprintf(fp,q的精度:%lf\n,mi[3]);fprintf(fp,w的精度:%lf\n,mi[4]);fprintf(fp,k的精度:%lf\n,mi[5]);}七、运算结果1.运行结果2.文件保存结果八、实习心得