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

《摄影测量学单像空间后方交会编程实习报告.doc》由会员分享，可在线阅读，更多相关《摄影测量学单像空间后方交会编程实习报告.doc（10页珍藏版）》请在装配图网上搜索。

摄影测量学 单像空间后方交会编程 实习报告 班 级： 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.99 36589.41 25273.32 2195.17 2 -53.40 82.21 37631.08 31324.51 728.69 3 -14.78 -76.63 39100.97 24934.98 2386.50 4 10.46 64.43 40426.54 30319.81 757.31 四、 实习原理 如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。 单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。 五、 实习流程 1. 获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。 2. 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。 3. 确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值： Zs0=H=m∙f Xs0=1ni-1nXti Ys0=1ni-1nYti φ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个未知数的中误差为： mi=Qiim0 当参加空间后方交会的控制点有n个时，则单位权中误差可按下式计算： m0=vv2n-6 流程图如下: 六、 程序代码 #include #include "Matrix.h" //矩阵运算头文件，来自网络，已上传 void main() { double Xs,Ys,Zs,q,w,k; //外方位元素 double x0,y0,f; //内方位元素 double x[4],y[4]; //影像坐标 double X0[4],Y0[4],Z0[4]; //地面坐标 double m; //比例尺 double a1,a2,a3,b1,b2,b3,c1,c2,c3; //旋转矩阵 int n=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 (int i=0;i<4;i++) { Xs += X0[i]; Ys += Y0[i]; } Xs /= 4; Ys /= 4; Zs=f*m; Matrix A(8, 6); Matrix L(8, 1); Matrix AT(6, 8); Matrix ATA(6, 6); Matrix ATA_(6, 6); Matrix Xx(6, 1); //迭代计算 do { //旋转矩阵R 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); //计算像点坐标 for (int i = 0; i < 4; i++) { double X = a1 * ( X0[i] - Xs) + b1 * (Y0[i] - Ys) + c1 * (Z0[i] - Zs); double Y = a2 * ( X0[i] - Xs) + b2 * (Y0[i] - Ys) + c2 * (Z0[i] - Zs); double Z = a3 * ( X0[i] - Xs) + b3 * (Y0[i] - Ys) + c3 * (Z0[i] - Zs); double xx = x[i]-x0; double yy = 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)_ATL AT = 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) && n<100); //求中误差 double m0,mi[6]; //单位权中误差的绝对值m0以及第i个未知数的中误差mi m0 = sqrt(((A*Xx - L).Transpose()*(A*Xx - L)).ToDouble() / (2 * 4 - 6)); Matrix Q(6,6); Q = (A.Transpose()*A).Inverse(); for (int i = 0; i < 6; 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=%.2lf Ys=%.2lf Zs=%.2lf\nq=%.5lf w=%.5lf k=%.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=%.2lf Ys=%.2lf Zs=%.2lf\nq=%.5lf w=%.5lf k=%.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. 文件保存结果 八、 实习心得