空间后方交会程序

1、. 实验目的：掌握摄影测量空间前方交会的原理，利用电脑编程语言实现空间前方交 会外方位元素的解算 。.仪器用具及数据文件：电脑windows xp系统，编程软件VISUAL C+6.0，地面控制点在摄 影测量坐标系中的坐标及其像点坐标文件shuju.txt 。:.实验内容：单张影像的空间前方交会：利用地面控制点数据及相应像点坐标 根据共线方程反求影像的外方位元素。数学模型：共线条件方程式：xy求解过程:f a1(XXs)b1(YYs)c1(ZZs )a3(XXs )b3(YYs)c3(ZZs )f a2(XXs)b2(YYs)c2(ZZs)a3(XXs)b3( YYs)c3(ZZs)1获取数据

2、。从航摄资料中查取平均航高与摄影机主距；获取控 制点的地面测量坐标并转换为地面摄影测量坐标。2量测控制点的像点坐标并做系统改正。3确定未知数的初始值。在竖直摄影且地面控制点大致分布均匀的情 况下，按如下方法确定初始值，即：0 X 0 Y 0ZXs，丫I, zS mf-nnn0 = 3= K=0式中；m为摄影比例尺分母；n为控制点个数。4用三个角元素的初始值，计算个方向余弦，组成旋转矩阵F。5逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面 坐标代入共线方程式，逐点计算像点坐标的近似值x y。6逐点计算误差方程式的系数和常数项，组成误差方程式。7计算法方程的系数矩阵 AtA和常数项AT

3、l，组成法方程式。8解法方程，求得外方位元素的改正数 dXs , dYs，dZs，d 0,d w,d k。9用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。第1页共8页10将求得的外方位元素改正数与规定的限差比拟，假设小于限差那么迭代结束。否那么用新的近似值重复49，直到满足要求为止。四实验结果 ：程序的源代码如下所示：#include<stdio.h>#include<stdlib.h> #include<malloc.h>#include<math.h> #include<conio.h> #define N 4

4、void 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+)

5、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&l

6、t;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

7、(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+)

8、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,

9、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"

10、;,&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);/

11、 输入 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)

12、;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-Y s0)+c2*(Zi-Zs0); / 计算像点坐标近似值y0i=-f*(a1*(Xi-Xs0)+b1*(Yi-Ys0)+c1*(Zi-Zs0)/(a2*(Xi-

13、Xs0)+b2*(Yi-Y s0)+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

14、=(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 ma

15、trix: ln");/ printmatrix(l,2*N,1); turn(A,AT,2*N,6);/ printf("output matrix: ATn");/ printmatrix(AT,6,2*N);mulAB(AT,A,A TA,6,2*N,2*N,6); / 间接平差法计算外方位元素 ,中间计算的矩阵可以根据需要 / printf("output matrix: ATAn");/选择性的输出，这里将其屏蔽，为了在打印输出结果的时候/ printmatrix(ATA,6,6); / 节约资源ATA_=inv(ATA,6);/ p

16、rintf("output matrix: ATA_n");f¥s*27483.39112, Zs =7593.2&e74rt0_ 00383,.0B04S ,<8ss sm号 key ta co nt inueH86.150000 r53-400600 kl4_7800B0-76.63000037631.9808003910H-?7000025273.320000 31324-51000024934.980000 30319-8100B02195 170000728.698S0B757.31G00B请输入摄影机主夢噂誌匕例尺分晞f ,n： 153-

17、24,50000Xs-397»9.Vs =27482-39112,26874,t0_002ea,w=-0.0Ha4«,k0.0G71£外方他元素为】/ prin tmatrix(ATA_,6,6);mulAB(AT,l,A TI,6,2*N,2*N,1);/ printf("outpit matrinx: ATIn);/ prin tmatrix(ATl,6,1);mulAB(ATA_,ATI,V ,6,6,6,1);/prin tf("output matrix: Vn );/prin tmatrix(V,6,1);XsO+=VO;YsO+=

18、V1;ZsO+=V2;t+=V3;w+=V4;k+=V5;printf("第d次计算的外方位元素为：n,+j);prin tf("Xs=%.5lf,Ys=%.5lf,Zs=%.5lf,t=%.5lf,w=%.5lf,k=%.5lfn",Xs0,Ys0,Zs0,t,w,k); printf("n外方位元素为：n");prin tf("Xs=%.5lf,Ys=%.5lf,Zs=%.5lf,t=%.5lf,w=%.5lf,k=%.5lfn",Xs0,Ys0,Zs0,t,w,k); fclose(fp);程序运行的结果为：飞八陳履V：语言源程序*貫世测量空阎厉方交会程序设计AMbu臥霞路测量仝阎前方交会 日回 旦|Xs=39787.81735,Vs=27612.87663,Zs=7545.60700,t=0亠炉二.側債