《空间后方交会编程实习报告(总5页)》由会员分享，可在线阅读，更多相关《空间后方交会编程实习报告(总5页)（5页珍藏版）》请在金锄头文库上搜索。
1、空间后方交会编程实习报告 一 实习目的用程序设计语言（Visual C+或者C语言）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。二 实习内容 利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。 三 实习数据 已知航摄仪的内方位元素：fk=153.24mm，x0=y0=0.0mm，摄影比例尺为1:500
2、00； 4个地面控制点的地面坐标及其对应像点的像片坐标：点号 像片坐标(mm) 地面点坐标(m) x y X Y Z 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个外方位元素，就能确定被摄物体与航摄影像的关系。因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。可采
3、取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，t，w，k。五 实习流程（1）获取已知数据。从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，f；获取控制点的空间坐标Xt
4、，Yt，Zt。（2）量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。（3）确定未知数的初始值。单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：Z0s=H=m*f+Zi/4；X0s=xi/n;Y0s=Yi/n;t=0;式中：m为摄影比例尺分母；(4)计算旋转矩阵R。利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。（5）逐点计算像点坐标的近似值。利用未知数的近似值按共线条件方程计算控制点像点坐标的近似值（x）、（y）。（6）按下式逐点计算误差方程式的系数和常数
5、项，组成误差方程。 (7)计算法方程的系数矩阵ATA与常数项ATL，组成法方程；（8）解求外方位元素。根据法方程，解求外方位元素的改正数，并与相应的近似值求和，得到外方位元素新的近似值。（9）检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较，通常对t、Xs、Ys、Zs的改正数t，Xs，Ys，Zs给予限差，当改正数小于限差时，迭代结束。否则用新的近似值重复（4）（8）步骤计算，直到满足要求为止。（10）空间后方交会的精度估计：按上述方法所求得的影像外方位元素的精度可以通过法方程式中未知数的系数矩阵的逆阵（ATA）-1来解求，此时视像点坐标为等精度不相关观测值。因为ATA）-1中第i
6、个主对角线上的元素Qii就是法方程式中第i个未知数的权倒数，若单位权中误差为m0，则第i个未知数的中误差为：mi=0当参加空间后方交会的控制点有n个时，则单位权中误差可按下式计算：六 主要代码与详解void R(double t,double w,double k,double *a,double *b,double *c) a0=cos(t)*cos(k)-sin(t)*sin(w)*sin(k); a1=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k); a2=-sin(t)*cos(w); b0=cos(w)*sin(k); b1=cos(w)*cos(k); b2
7、=-sin(w); c0=sin(t)*cos(k)+cos(t)*sin(w)*sin(k); c1=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k); c2=cos(t)*cos(w);/子函数计算旋转矩阵R。利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。void main() int i,m,num;double t,w,k,Xs,Ys,Zs,f; /六个外方位元素与焦距f double xN=-86.15,-53.40,-14.78,10.46,yN=-68.99,82.21,-76.63,64.43;/4个
8、像点坐标double XN=36589.41,37631.08,39100.97,40426.54,YN=25273.32,31324.51,24934.98,30319.81,ZN=2195.17,728.69,2386.50,757.31;/四个控制点的空间坐标m=50000;f=153.24;/ 影像比例尺1/m 与焦距f/ 以上主要是已知数据的定义double V6=0,a3,b3,c3,XoN,YoN,ZoN,A6*M,B6*M,lM,C36,D6,E8;for(i=0;iN;i+) xi=xi/1000.0; yi=yi/1000.0;/像点坐标单位换算 / 以下为确定未知数的初始
9、值。 Xs=(X0+X1+X2+X3)/N; Ys=(Y0+Y1+Y2+Y3)/N; Zs=m*f+(Z0+Z1+Z2+Z3)/N;/ Xs，Ys，Zs 为摄站点的空间坐标初始值 f=f/1000.0; t=w=k=0.0; /角元素的近似值 for(num=1;num10;num+) /num控制循环迭代次数R(t,w,k,a,b,c);/ 计算旋转矩阵R。利用角元素的近似值按下式计算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。 for(i=0;iN;i+)/下面的循环用来逐点计算误差方程式的系数和常数项，组成误差方程。 /以下是用共线条件方程计算控制点像点坐标
10、的近似值（x）、（y）。 Xoi=-f*(a0*(Xi-Xs)+b0*(Yi-Ys)+c0*(Zi-Zs)/(a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs); Yoi=-f*(a1*(Xi-Xs)+b1*(Yi-Ys)+c1*(Zi-Zs)/(a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs);Zoi=a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs);/Zo便于后面计算 A12*i+0=(a0*f+a2*xi)/Zoi; A12*i+1=(b0*f+b2*xi)/Zoi; A12*i+2=(c0*f+c2*xi)/Zoi; A12*i+3=yi*
11、sin(w)-(xi*(xi*cos(k)-yi*sin(k)/f+f*cos(k)*cos(w); A12*i+4=-f*sin(k)-xi*(xi*sin(k)+yi*cos(k)/f; A12*i+5=yi; A12*i+6=(a1*f+a2*yi)/Zoi; A12*i+7=(b1*f+b2*yi)/Zoi; A12*i+8=(c1*f+c2*yi)/Zoi; A12*i+9=-xi*sin(w)-(yi*(xi*cos(k)-yi*sin(k)/f-f*sin(k)*cos(w); A12*i+10=-f*cos(k)-yi*(xi*sin(k)+yi*cos(k)/f; A12*i
12、+11=-xi; l2*i=xi-Xoi; l2*i+1=yi-Yoi;/计算l zhuanzhi(A,B,8,6); /转置A xiangchen(B,A,C,6,8,6);/求ATA xiangchen(B,l,D,6,8,1); /求常数项ATL qiuni(C,6);/求ATA的逆 xiangchen(C,D,V,6,6,1);/求改正数V Xs+=V0; Ys+=V1; Zs+=V2; t+=V3; w+=V4; k+=V5; /结果改正 if(fabs(V0)=0.00001)&(fabs(V1)=0.00001)&(fabs(V2)=0.00001)&(fabs(V3)=0.00
13、001)&(fabs(V4)=0.00001)&(fabs(V5)=0.00001) break; /限差控制。检查计算是否收敛。将所求得的外方位元素的改正数与规定的限差比较。当改正数小于限差时，迭代结束。否则用新的近似值重复计算，直到满足要求为止 。 /以下是计算单位权中误差double s=0,m0,v8; xiangchen(A,V,E,8,6,1); /计算AX for(i=0;i8;i+) vi=Ei-li;/计算AX-L s+=vi*vi;/计算vv m0=sqrt(s/2);/ 单位权中误差按下式计算：，这里n=4. R(t,w,k,a,b,c);/输出结果 printf(像主点的空间坐标为：n); printf(Xs=%.2fn,Xs); printf(Ys=%.2fn,Ys); printf(Zs=%.2fn,Zs); printf(旋转矩阵R为:n); for(i=0;i3;i+) printf(
