摄影测量学单像空间后方交会编程实习报告材料.doc
word摄影测量学单像空间后方交会编程实习报告一、 实习目的通过对提供的数据进展计算，输出像片的外方位元素并评定精度。深入理解单像空间后方交会的思想，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过尝试编程实现加强编程处理问题的能力和对实习内容的理解，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。了解摄影测量平差的根本过程，掌握空间后方交会的定义和实现算法。二、 实习内容根据学习的单像空间后方交会的知识，用程序设计语言C+或C语言编写一个完整的单像空间后方交会程序，通过对提供的数据进展计算，输出像片的外方位元素并评定精度。三、 实习数据航摄仪的内方位元素：fk=153.24mm，x0=y0=0，摄影比例尺为1:15000；4个地面控制点的地面坐标与其对应像点的像片坐标： 像片坐标地面点坐标x(mm)y(mm)X(m)Y(m)Z(m)1234四、 实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采取的方法有：利用雷达、全球定位系统GPS、惯性导航系统INS以与星相摄影机来获取影像的外方位元素；也可以利用影像覆盖X围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。 单像空间后方交会的根本思想是：以单幅影像为根底，从该影像所覆盖地面X围内假如干控制点的地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,。五、 实习流程1. 获取数据。从摄影资料中查取影像比例尺1/m，平均摄影距离航空摄影的航高、内方位元素x0，y0，；获取控制点的空间坐标Xt，Yt，Zt。 2. 量测控制点的像点坐标并进展必要的影像坐标系统误差改正，得到像点坐标。3. 确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：式中：m为摄影比例尺分母,n为控制点个数； 4. 计算旋转矩阵R。利用角元素的近似值按公式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。 5. 逐点计算像点坐标的近似值。利用未知数的近似值按共线条件方程计算控制点像点坐标的近似值(x),(y)。 6. 逐点计算误差方程式的系数和常数项，组成误差方程。7. 计算法方程的系数矩阵ATA与常数项ATL，组成法方程。 8. 解求外方位元素。根据法方程，解求外方位元素的改正数，并与相应的近似值求和，得到外方位元素新的近似值。9. 检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比拟，通常对,的改正数给予限差，当改正数小于限差时，迭代完毕。否如此用新的近似值重复48步骤计算，直到满足要求为止。 10. 空间后方交会的精度估计： 按上述方法所求得的影像外方位元素的精度可以通过法方程式中未知数的系数矩阵的逆阵ATA-1来解求，此时视像点坐标为等精度不相关观测值。因为(ATA-1中第i个主对角线上的元素Qii就是法方程式中第i个未知数的权倒数，假如单位权中误差为m0，如此第i个未知数的中误差为：当参加空间后方交会的控制点有n个时，如此单位权中误差可按下式计算：流程图如下:六、 程序代码#include <stdio.h>#include "Matrix.h" /矩阵运算头文件，来自网络，已上传void main()double Xs,Ys,Zs,q,w,k;/外方位元素double x0,y0,f;/内方位元素double x4,y4;/影像坐标double X04,Y04,Z04; /地面坐标double m;/比例尺double a1,a2,a3,b1,b2,b3,c1,c2,c3; /旋转矩阵int n=0;/迭代次数/地面坐标X00=36589.41;X01=37631.08;X02=39100.97;X03=40426.54;Y00=25273.32;Y01=31324.51;Y02=24934.98;Y03=30319.81;Z00=2195.17;Z01=728.69;Z02=2386.50;Z03=757.31;/影像坐标x0=-0.08615;x1=-0.05340;x2=-0.01478;x3=0.01046;y0=-0.06899;y1=0.08221;y2=-0.07663;y3=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 += X0i;Ys += Y0i;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/旋转矩阵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 (int i = 0; i < 4; i+)double X = a1 * ( X0i - Xs) + b1 * (Y0i - Ys) + c1 * (Z0i - Zs);double Y = a2 * ( X0i - Xs) + b2 * (Y0i - Ys) + c2 * (Z0i - Zs);double Z = a3 * ( X0i - Xs) + b3 * (Y0i - Ys) + c3 * (Z0i - Zs);double xx = xi-x0;double yy = yi-y0;A2 * i0 = (a1*f + a3*xx) / Z;A2 * i1 = (b1*f + b3*xx) / Z;A2 * i2 = (c1*f + c3*xx) / Z;A2 * i3 = yy*sin(w) - (xx*(xx*cos(k) - yy*sin(k) / f + f*cos(k)*cos(w);A2 * i4 = -f*sin(k) - xx*(xx*sin(k) + yy*cos(k) / f;A2 * i5 = yy;A2 * i + 10 = (a2*f + a3*yy) / Z;A2 * i + 11 = (b2*f + b3*yy) / Z;A2 * i + 12 = (c2*f + c3*yy) / Z;A2 * i + 13 = -xx*sin(w) - (yy*(xx*cos(k) - yy*sin(k) / f - f*sin(k)*cos(w);A2 * i + 14 = -f*cos(k) - yy*(xx*sin(k) + yy*cos(k) / f;A2 * i + 15 = -xx;L2 * i0 = xi - (-f*X / Z);L2 * i + 10 =yi - (-f*Y / Z);/法方程的解X=(ATA)_ATLAT = A.Transpose(); /转置ATA = AT * A;ATA_ = ATA.Inverse(); /求逆Xx = (ATA_ * AT) * L; /法方程解Xs += Xx00;Ys += Xx10;Zs += Xx20;q += Xx30;w += Xx40;k += Xx50;n+;/printf("%dn",n);/printf("%lf %lf %lf %.5lf %.5lf %.5lfn", Xs, Ys, Zs, q, w, k);while (abs(Xx30) > 0.000001 | abs(Xx40) > 0.000001 | abs(Xx50) > 0.000001) && n<100);/求中误差double m0,mi6; /单位权中误差的绝对值m0以与第i个未知数的中误差mim0 = 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+)mii = (sqrt(Qii)*m0;/printf("%lfn",m0);/输出printf("迭代次数:%dn",n);printf("旋转矩阵R:n");printf("%.5lf %.5lf %.5lfn",a1,a2,a3);printf("%.5lf %.5lf %.5lfn",b1,b2,b3);printf("%.5lf %.5lf %.5lfn",c1,c2,c3);printf("n外方位元素解:nWs=%.2lf Ys=%.2lf Zs=%.2lfnq=%.5lf w=%.5lf k=%.5lfnn", Xs, Ys, Zs, q, w, k);printf("单位权中误差的绝对值:%lfmn",m0);printf("Xs的精度:%lfmn",mi0);printf("Ys的精度:%lfmn",mi1);printf("Zs的精度:%lfmn",mi2);printf("q的精度:%lfn",mi3);printf("w的精度:%lfn",mi4);printf("k的精度:%lfn",mi5);/保存到文件，结果文件默认保存在程序目录FILE* fp;fp=fopen("result.txt","w");/fprintf(fp,"迭代次数:%dn",n);fprintf(fp,"旋转矩阵R:n");fprintf(fp,"%.5lf %.5lf %.5lfn",a1,a2,a3);fprintf(fp,"%.5lf %.5lf %.5lfn",b1,b2,b3);fprintf(fp,"%.5lf %.5lf %.5lfn",c1,c2,c3);fprintf(fp,"n外方位元素解:nWs=%.2lf Ys=%.2lf Zs=%.2lfnq=%.5lf w=%.5lf k=%.5lfnn", Xs, Ys, Zs, q, w, k);fprintf(fp,"单位权中误差的绝对值:%lfmn",m0);fprintf(fp,"Xs的精度:%lfmn",mi0);fprintf(fp,"Ys的精度:%lfmn",mi1);fprintf(fp,"Zs的精度:%lfmn",mi2);fprintf(fp,"q的精度:%lfn",mi3);fprintf(fp,"w的精度:%lfn",mi4);fprintf(fp,"k的精度:%lfn",mi5);七、 运算结果1. 运行结果2. 文件保存结果13 / 13
