空间后方交会程序.doc
一 实验目的： 掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交会外方位元素的解算。二 仪器用具与数据文件： 计算机windows xp 系统，编程软件VISUAL C+6.0，地面控制点在摄影测量坐标系中的坐标与其像点坐标文件shuju.txt。三 实验内容： 单影像的空间后方交会：利用地面控制点数据与相应像点坐标根据共线方程反求影像的外方位元素。 数学模型：共线条件方程式：求解过程：1获取数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影测量坐标。2量测控制点的像点坐标并做系统改正。3确定未知数的初始值。在竖直摄影且地面控制点大致分布均匀的情况下，按如下方法确定初始值，即：, =0 式中；m为摄影比例尺分母；n为控制点个数。 4用三个角元素的初始值，计算个方向余弦，组成旋转矩阵R。 5逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值x、y。 6逐点计算误差方程式的系数和常数项，组成误差方程式。7计算法方程的系数矩阵和常数项，组成法方程式。 8解法方程，求得外方位元素的改正数,，d,d,d。9用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。10将求得的外方位元素改正数与规定的限差比拟，假如小于限差如此迭代完毕。否如此用新的近似值重复49，直到满足要求为止。四实验结果：程序的源代码如下所示：#include<stdio.h>#include<stdlib.h>#include<malloc.h>#include<math.h>#include<conio.h>#define N 4void turn(double *A,double A2,int m,int n) /计算矩阵的转置int i,j; for(i=0;i<m;i+)for(j=0;j<n;j+)A2j*m+i=Ai*n+j;void mulAB(double *A,double *B,double *C,int am,int an,int bm,int bn) /计算两矩阵相乘int i,j,l,u;if(an!=bm)printf("error!cannot do the multiplication.n");return;for(i=0;i<am;i+)for(j=0;j<bn;j+)u=i*bn+j;Cu=0.0;for(l=0;l<an;l+)Cu+=Ai*an+l*Bl*bn+j;return;double *inv(double *a,int n) /计算矩阵的逆，本程序的难点 /采用高斯-约旦-全选主元法int *is,*js,i,j,k,l,u,v; double d,p; is=(int*)malloc(n*sizeof(int); js=(int*)malloc(n*sizeof(int); for (k=0; k<=n-1; k+) d=0.0; for (i=k;i<n;i+) for (j=k;j<n;j+) l=i*n+j; p=fabs(al); if (p>d) d=p; isk=i; jsk=j; if (d+1.0=1.0) free(is); free(js); printf("error not invn"); return NULL; if (isk!=k) for (j=0;j<n;j+) u=k*n+j; v=isk*n+j; p=au; au=av; av=p; if (jsk!=k) for (i=0;i<n;i+) u=i*n+k; v=i*n+jsk; p=au; au=av; av=p; l=k*n+k; al=1.0/al; for (j=0;j<n;j+) if (j!=k) u=k*n+j; au=au*al; for (i=0;i<n;i+) if (i!=k) for (j=0;j<n;j+) if (j!=k) u=i*n+j; au=au-ai*n+k*ak*n+j; for (i=0;i<n;i+) if (i!=k) u=i*n+k; au=-au*al; for (k=n-1;k>=0;k-) if (jsk!=k) for (j=0;j<=n-1;j+) u=k*n+j; v=jsk*n+j; p=au; au=av; av=p; if (isk!=k) for (i=0;i<n;i+) u=i*n+k; v=i*n+isk; p=au; au=av; av=p; free(is);free(js); return a;void printmatrix(double M,int m,int n) /矩阵的输出int i,j;for(i=0;i<m;i+)for(j=0;j<n;j+)printf("%.5f",Mi*n+j);printf("n");printf("n");main() /主函数，空间后方交会的计算FILE *fp; /定义一个文件指针fpint m,i,j=0;double f,t,w,k,S1=0.0,S2=0.0,S3=0.0,xN,yN,x0N,y0N,XN,YN,ZN,Xs0,Ys0,Zs0;double a3,b3,c3,A2*N*6,AT6*2*N,ATA6*6,*ATA_=NULL,l2*N,ATl6,V6; if(fp=fopen("e:shuju.txt","r")=NULL) / /并判断文件是否打开成功printf("nerror on open shuju.txtn");getch();exit(1); for(i=0;i<N;i+)fscanf(fp,"%lf%lf%lf%lf%lf",&xi,&yi,&Xi,&Yi,&Zi); /将文件中的数据赋给主函数定义的变量printf("原始数据:n");printf("xttyttXttYttZttn"); /输出文件中的原始数据for(i=0;i<N;i+)printf("%lf %lf %lf %lf %lfn",xi,yi,Xi,Yi,Zi);printf("n请输入摄影机主距和摄影比例尺分母；f,m:");scanf("%lf,%lf",&f,&m); /输入f,mf=f/1000.0;for(i=0;i<N;i+)xi/=1000.0;yi/=1000.0;S1+=Xi;S2+=Yi;S3+=Zi; Xs0=S1/N;Ys0=S2/N; /计算外方位元素的初始值Zs0=m*f+S3/N;t=0.0;w=0.0;k=0.0; while(j<3) printf("-第%d次计算-n",j+1);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=-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);for(i=0;i<N;i+)x0i=-f*(a0*(Xi-Xs0)+b0*(Yi-Ys0)+c0*(Zi-Zs0)/(a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0); /计算像点坐标近似值y0i=-f*(a1*(Xi-Xs0)+b1*(Yi-Ys0)+c1*(Zi-Zs0)/(a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0);Gi=a2*(Xi-Xs0)+b2*(Yi-Ys0)+c2*(Zi-Zs0);for(i=0;i<N;i+)Ai*12+0=(a0*f+a2*xi)/Gi;Ai*12+1=(b0*f+b2*xi)/Gi;Ai*12+2=(c0*f+c2*xi)/Gi;Ai*12+3=yi*sin(w)-(xi*(xi*cos(k)-yi*sin(k)/f+f*cos(k)*cos(w); /计算误差方程的系数阵以与lAi*12+4=-f*sin(k)-xi*(xi*sin(k)+yi*cos(k)/f;Ai*12+5=yi;Ai*12+6=(a1*f+a2*yi)/Gi;Ai*12+7=(b1*f+b2*yi)/Gi;Ai*12+8=(c1*f+c2*yi)/Gi;Ai*12+9=-xi*sin(w)-(yi*(xi*cos(k)-yi*sin(k)/f-f*sin(k)*cos(w);Ai*12+10=-f*cos(k)-yi*(xi*sin(k)+yi*cos(k)/f;Ai*12+11=-xi;li*2+0=xi-x0i;li*2+1=yi-y0i;/printf("output matrix: An");/printmatrix(A,2*N,6);/printf("output matrix: ln");/ printmatrix(l,2*N,1);turn(A,AT,2*N,6);/printf("output matrix: ATn");/printmatrix(AT,6,2*N);mulAB(AT,A,ATA,6,2*N,2*N,6); /间接平差法计算外方位元素,中间计算的矩阵可以根据需要/printf("output matrix: ATAn"); /选择性的输出，这里将其屏蔽，为了在打印输出结果的时候/printmatrix(ATA,6,6); /节约资源ATA_=inv(ATA,6);/printf("output matrix: ATA_n");/printmatrix(ATA_,6,6);mulAB(AT,l,ATl,6,2*N,2*N,1);/printf("outpit matrinx: ATln");/printmatrix(ATl,6,1);mulAB(ATA_,ATl,V,6,6,6,1);/printf("output matrix: Vn");/printmatrix(V,6,1);Xs0+=V0;Ys0+=V1;Zs0+=V2;t+=V3;w+=V4;k+=V5;printf("第%d次计算的外方位元素为：n",+j); printf("Xs=%.5lf,Ys=%.5lf,Zs=%.5lf,t=%.5lf,w=%.5lf,k=%.5lfn",Xs0,Ys0,Zs0,t,w,k); printf("n外方位元素为：n"); printf("Xs=%.5lf,Ys=%.5lf,Zs=%.5lf,t=%.5lf,w=%.5lf,k=%.5lfn",Xs0,Ys0,Zs0,t,w,k); fclose(fp);程序运行的结果为：五实验总结： 通过这次的实验我学到了很多的东西，通过编程加深了对摄影测量空间后方交会相关知识的理解。在教师的严格要求下，我翻阅了很多的C语言的书籍，看了很多的算法，这个程序最难的地方是矩阵的求逆，以与如何将文本文档中的数据赋给程序中的数组变量，这些在以前的学习中都没有学到，感谢昭华教师的严格要求。
