前方后方空间交会实验报告

1、.中南大学本科生课程设计(实践)任务书、设计报告 (摄影测量与遥感概论)题 目空间后方-前方交会 学生姓名指导教师邹峥嵘学 院地球科学与信息物理学院专业班级测绘0902班学生学号一、 实验目的通过对数字影像空间后交前交的程序设计实验，要求我们进一步理解和掌握影像外方位元素的有关理论、原理和方法。利用计算机程序设计语言编写摄影测量空间交会软件进行快速确定影像的外方位元素及其精度，然后通过求得的外方位元素求解未知点的地面摄影测量坐标，达到通过摄影测量量测地面地理数据的目的。二、 实验要求Ø 用C、VB或者Matlab编写空间后方交会-前方交会计算机程序。Ø 提交实验报告：程序框

2、图，程序源代码、计算结果及体会。Ø 计算结果：地面点坐标、外方位元素及精度。Ø 完成时间：2011年11月17日。三、 实验数据点号 左片 右片 地面摄影测量坐标 x y x y X Y Z GCP1 16.012 79.963 -73.93 78.706 5083.205 5852.099 527.925 GCP2 88.56 81.134 -5.252 78.184 5780.02 5906.365 571.549 GCP3 13.362 -79.37 -79.122 -78.879 5210.879 4258.446 461.81 GCP4 82.24 -80.027

3、 -9.887 -80.089 5909.264 4314.283 455.484 1 51.758 80.555 -39.953 78.463 2 14.618 -0.231 -76.006 0.036 3 49.88 -0.782 -42.201 -1.022 4 86.14 -1.346 -7.706 -2.112 5 48.035 -79.962 -44.438 -79.736 f=150.000mm，x0=0，y0=0四、 实验思路Ø 利用后方交会得出两张像片各自的外方位元素1) 获取已知数据：从摄影资料中插曲像片比例尺、平均航高、内方位元素以及控制点的地面摄影测量坐标及对

4、应的像点坐标。2) 确定未知数的初始值：在竖直摄影的情况下，胶原素的初始值为0，线元素其中Zs=m*f+，Xs=，Ys=。3) 计算旋转矩阵R。4) 逐点计算像点坐标的近似值：利用共线方程。5) 组成误差方程并法化。6) 解求外方位元素。7) 检查计算是否收敛。Ø 利用解求出的外方位元素进行前方交会1) 用各自像片的角元素计算出左右像片的旋转矩阵R1和R2。2) 根据左右像片的外方位元素计算摄影基线分量Bx,By，Bz。3) 逐点计算像点的空间辅助坐标。4) 计算投影系数。5) 计算未知点的地面摄影测量坐标。6) 重复以上步骤完成所有点的地面坐标的计算。五、 实验过程Ø 程

5、序流程框图后方交会函数确定已知数据比例尺m确定各外方位元素初始值计算旋转矩阵逐点计算像点坐标近似值不满足限差则重复计算逐点计算误差方程系数项，组成误差系数矩阵A利用矩阵运算求解外方位元素检查是否满足限差若满足则输出外方位元素将整个过程作为一个函数继续进行右片的外方位元素求解求解各外方位元素精度空间前方交会利用已求得的角元素计算2张像片各自的旋转矩阵利用已求得的线元素Xs1，Ys1，Zs1，p0，,w01，k01；Xs2，Ys2，Zs2，p02，w02，k02，计算基线分量：Bx=Xs2-Xs1；By=Ys2-Ys1；Bz=Zs2-Zs1；输入像片坐标，利用旋转矩阵求解想空间辅助坐标计算点投影系

6、数：N1=(Bx*Z2-Bz*X2)/(X1*Z2-X2*Z1)；N2=(Bx*Z1-Bz*X1)/(X1*Z2-X2*Z1)；计算地面摄影测量坐标Xt=(N1*X1+Xs1)+(N2*X2+Xs2)/2；Yt=(N1*Y1+Ys1)+(N2*Y2+Ys2)/2；Zt=(N1*Z1+Zs1)+(N2*Z2+Zs2) /2 ；结束程序Ø 程序中的主要函数设计子函数（矩阵求积multiply，计算函数Resection，矩阵转置transpose，矩阵求逆inMerse1，输出函数shuchu，左片的外方位元素求解函数zuobian。右片的外方位元素求解函数youbian。）Ø

7、 程序源代码#include "stdio.h"#include "math.h"double Xs1,Xs2,Ys1,Ys2,Zs1,Zs2,p01,p02,w01,w02,k01,k02;/求矩阵a的转置矩阵b，a为m行、n列void transpose(double *a, double *b, int m, int n);/矩阵a乘以矩阵b，结果存储在c中，a为m×n大小，b为n×l大小void multiply(double *a, double *b, double *c, int m, int n, int l);/求矩

8、阵a的逆int inMerse1(double *a, int n);/输出m行、n列的矩阵avoid shuchu(double *a, int m, int n);/计算并输出左片的外方位元素void zuobian();/计算并输出右片的外方位元素void youbian();void zuobian() FILE *fp = NULL;FILE *fp1 = NULL;if(fp=fopen("F:image.txt","r") = NULL)printf("Open file error!");return;if(fp1=f

9、open("F:ground.txt","r") = NULL)printf("Open file error!");return;/像点坐标和地面点坐标double imagecontrol42=0.0;double groundcontrol43=0.0;/摄影比例尺分母double m = 9943;double f=0.15;long i,j,k; for(i=0; i<4; i+)for(j=0; j<2; j+)fscanf(fp, "%lf", &imagecontrolij);i

10、magecontrolij /= 1000.0;for(k=0; k<3; k+)fscanf(fp1, "%lf", &groundcontrolik);fclose(fp);fclose(fp1); /计算外方位元素初始值for( i=0;i<4;i+)Xs1+=groundcontroli0;Ys1+=groundcontroli1;Zs1+=groundcontroli2;Xs1/=4.0;Ys1/=4.0; Zs1/=4.0;Zs1+=m*f;double R33=0.0;double L3=0.0,L1=0.0,L2=0.0;double L

11、81=0.0,x=0.0,y=0.0;double A86=0.0,AT68=0.0,ATA66=0.0,B68=0.0;double V61=0.0;int n=0;do/计算旋转矩阵R00=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);R01=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R02=(-1)*sin(p01)*cos(w01);R10=cos(w01)*sin(k01);R11=cos(w01)*cos(k01);R12=(-1)*sin(w01);R20=sin(p01)*co

12、s(k01)+cos(p01)*sin(w01)*sin(k01);R21=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R22=cos(p01)*cos(w01);for(i=0,j=0;j<4;i+=2,j+)/计算像点坐标的近似值L1=R00*(groundcontrolj0-Xs1)+R10*(groundcontrolj1-Ys1)+R20*(groundcontrolj2-Zs1);L2=R01*(groundcontrolj0-Xs1)+R11*(groundcontrolj1-Ys1)+R21*(groundcontr

13、olj2-Zs1);L3=R02*(groundcontrolj0-Xs1)+R12*(groundcontrolj1-Ys1)+R22*(groundcontrolj2-Zs1);x=(-1)*f*L1/L3;y=(-1)*f*L2/L3; /计算常数项L2*j0=imagecontrolj0-x;L2*j+10=imagecontrolj1-y;/计算系数矩阵Ai0=(R00*f+R02*imagecontrolj0)/L3;Ai1=(R10*f+R12*imagecontrolj0)/L3;Ai2=(R20*f+R22*imagecontrolj0)/L3;Ai3=imagecontro

14、lj1*sin(w01)-(imagecontrolj0/f)*(imagecontrolj0*cos(k01)-imagecontrolj1*sin(k01)+f*cos(k01)*cos(w01);Ai4=(-1)*f*sin(k01)-(imagecontrolj0/f)*(imagecontrolj0*sin(k01)+imagecontrolj1*cos(k01);Ai5=imagecontrolj1;Ai+10=(R01*f+R02*imagecontrolj1)/L3;Ai+11=(R11*f+R12*imagecontrolj1)/L3;Ai+12=(R21*f+R22*ima

15、gecontrolj1)/L3;Ai+13=(-1)*imagecontrolj0*sin(w01)-(imagecontrolj1/f)*(imagecontrolj0*cos(k01)-imagecontrolj1*sin(k01)-f*sin(k01)*cos(w01);Ai+14=(-1)*f*cos(k01)-(imagecontrolj1/f)*(imagecontrolj0*sin(k01)+imagecontrolj1*cos(k01);Ai+15=(-1)*imagecontrolj0;transpose(&A00,&AT00,8,6);multiply(&a

16、mp;AT00,&A00,&ATA00,6,8,6);inMerse1(*ATA,6); multiply(&ATA00,&AT00,&B00,6,6,8);multiply(&B00,&L00,&V00,6,8,1); Xs1+=V00;Ys1+=V10;Zs1+=V20;p01+=V30;w01+=V40;k01+=V50;n+;while(fabs(V30)>=0.00001|fabs(V40)>=0.00001|fabs(V50)>=0.00001);R00=cos(p01)*cos(k01)-sin(p

17、01)*sin(w01)*sin(k01);R01=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R02=(-1)*sin(p01)*cos(w01);R10=cos(w01)*sin(k01);R11=cos(w01)*cos(k01);R12=(-1)*sin(w01);R20=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R21=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R22=cos(p01)*cos(w01);/进行未知数的精度评

18、定double AV81,X81,XT18,XTX11,mo,D66,mi6;multiply(&A00,&V00,&AV00,8,6,1);for(i=0;i<8;i+)Xi0=AVi0-Li0;transpose(&X00,&XT00,8,1);multiply(&XT00,&X00,&XTX00,1,8,1);mo=XTX00/2; for(i=0;i<6;i+)for(j=0;j<6;j+)Dij=mo*ATAij; for(i=0;i<6;i+)mii=sqrt(Dii);printf("

19、;左片结果为：nn");printf("旋转矩阵R为：n");shuchu(&R00,3,3);printf("外方为元素为：n");printf("Xs1=%lfn",Xs1);printf("Ys1=%lfn",Ys1);printf("Zs1=%lfn",Zs1);printf("p01=%lfn",p01);printf("w01=%lfn",w01);printf("k01=%lfn",k01);printf

20、("各外方位元素精度为：n"); printf("m1=%lfnm2=%lfnm3=%lfnm4=%lfnm5=%lfnm6=%lfn",mi0,mi1,mi2,mi3,mi4,mi5);printf("迭代次数为：%dnnnn",n);fclose(fp);void youbian()FILE *fp = NULL;FILE *fp1 = NULL;if(fp=fopen("F:image2.txt","r") = NULL)printf("Open file error!"

21、;);return;if(fp1=fopen("F:ground.txt","r") = NULL)printf("Open file error!");return;/像点坐标和地面点坐标double imagecontrol42=0.0;double groundcontrol43=0.0;/摄影比例尺分母double m = 10177;double f=0.15;long i,j,k; for(i=0; i<4; i+)for(j=0; j<2; j+)fscanf(fp, "%lf", &am

22、p;imagecontrolij);imagecontrolij /= 1000.0;for(k=0; k<3; k+)fscanf(fp1, "%lf", &groundcontrolik);fclose(fp);fclose(fp1); /计算外方位元素初始值for( i=0;i<4;i+)Xs2+=groundcontroli0;Ys2+=groundcontroli1;Zs2+=groundcontroli2;Xs2/=4.0;Ys2/=4.0; Zs2/=4.0;Zs2+=m*f;double R33=0.0;double L3=0.0,L1=

23、0.0,L2=0.0;double L81=0.0,x=0.0,y=0.0;double A86=0.0,AT68=0.0,ATA66=0.0,B68=0.0;double V61=0.0;int n=0;doR00=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R01=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R02=(-1)*sin(p02)*cos(w02);R10=cos(w02)*sin(k02);R11=cos(w02)*cos(k02);R12=(-1)*sin(w02);R20

24、=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R21=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R22=cos(p02)*cos(w02);for(i=0,j=0;j<4;i+=2,j+)/计算像点坐标的近似值L1=R00*(groundcontrolj0-Xs2)+R10*(groundcontrolj1-Ys2)+R20*(groundcontrolj2-Zs2);L2=R01*(groundcontrolj0-Xs2)+R11*(groundcontrolj1-Ys2)+R21*

25、(groundcontrolj2-Zs2);L3=R02*(groundcontrolj0-Xs2)+R12*(groundcontrolj1-Ys2)+R22*(groundcontrolj2-Zs2);x=(-1)*f*L1/L3;y=(-1)*f*L2/L3; /计算常数项矩阵L2*j0=imagecontrolj0-x;L2*j+10=imagecontrolj1-y;/计算系数矩阵Ai0=(R00*f+R02*imagecontrolj0)/L3;Ai1=(R10*f+R12*imagecontrolj0)/L3;Ai2=(R20*f+R22*imagecontrolj0)/L3;A

26、i3=imagecontrolj1*sin(w02)-(imagecontrolj0/f)*(imagecontrolj0*cos(k02)-imagecontrolj1*sin(k02)+f*cos(k02)*cos(w02);Ai4=(-1)*f*sin(k02)-(imagecontrolj0/f)*(imagecontrolj0*sin(k02)+imagecontrolj1*cos(k02);Ai5=imagecontrolj1;Ai+10=(R01*f+R02*imagecontrolj1)/L3;Ai+11=(R11*f+R12*imagecontrolj1)/L3;Ai+12=

27、(R21*f+R22*imagecontrolj1)/L3;Ai+13=(-1)*imagecontrolj0*sin(w02)-(imagecontrolj1/f)*(imagecontrolj0*cos(k02)-imagecontrolj1*sin(k02)-f*sin(k02)*cos(w02);Ai+14=(-1)*f*cos(k02)-(imagecontrolj1/f)*(imagecontrolj0*sin(k02)+imagecontrolj1*cos(k02);Ai+15=(-1)*imagecontrolj0;transpose(&A00,&AT00,8,

28、6);multiply(&AT00,&A00,&ATA00,6,8,6);inMerse1(*ATA,6); multiply(&ATA00,&AT00,&B00,6,6,8);multiply(&B00,&L00,&V00,6,8,1); Xs2+=V00;Ys2+=V10;Zs2+=V20;p02+=V30;w02+=V40;k02+=V50;n+;while(fabs(V30)>=0.00001|fabs(V40)>=0.00001|fabs(V50)>=0.00001);R00=cos(p02)*

29、cos(k02)-sin(p02)*sin(w02)*sin(k02);R01=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R02=(-1)*sin(p02)*cos(w02);R10=cos(w02)*sin(k02);R11=cos(w02)*cos(k02);R12=(-1)*sin(w02);R20=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R21=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R22=cos(p02)*cos(w

30、02);/进行未知数的精度评定double AV81,X81,XT18,XTX11,mo,D66,mi6;multiply(&A00,&V00,&AV00,8,6,1);for(i=0;i<8;i+)Xi0=AVi0-Li0;transpose(&X00,&XT00,8,1);multiply(&XT00,&X00,&XTX00,1,8,1);mo=XTX00/2; for(i=0;i<6;i+)for(j=0;j<6;j+)Dij=mo*ATAij; for(i=0;i<6;i+)mii=sqrt(Dii

31、);printf("右片结果为：nn");printf("旋转矩阵R为：n");shuchu(&R00,3,3);printf("外方为元素为：n");printf("Xs2=%lfn",Xs2);printf("Ys2=%lfn",Ys2);printf("Zs2=%lfn",Zs2);printf("p02=%lfn",p02);printf("w02=%lfn",w02);printf("k02=%lfn&quo

32、t;,k02);printf("各外方位元素精度为：n"); printf("m1=%lfnm2=%lfnm3=%lfnm4=%lfnm5=%lfnm6=%lfn",mi0,mi1,mi2,mi3,mi4,mi5);printf("迭代次数为：%dn",n);fclose(fp);/求矩阵a的转置矩阵b，a为m行、n列void transpose(double *a, double *b, int m, int n)int i,j;for(i=0;i<m;i+)for(j=0;j<n;j+)bj*m+i = ai*n+j;

33、/矩阵a乘以矩阵b，结果存储在c中，a为m×n大小，b为n×l大小void multiply(double *a, double *b, double *c, int m, int n, int l)int i,j,k;double t;for(i=0;i<m;i+)for(j=0;j<l;j+)t=0;for(k=0;k<n;k+)t += ai*n+k*bk*l+j;ci*l+j=t;/求矩阵a的逆int inMerse1(double *a, int n)int *is, *js, i, j, k, l, u, M;double d,p;is = n

34、ew intn;js = new intn;for(k=0; k<=n-1; k+)d=0.0;for(i=k; i<=n-1; i+)for(j=k; j<=n-1; j+)l=i*n+j; p=fabs(al);if (p>d) d=p; isk=i; jsk=j;if(d+1.0=1.0)delete is;delete js;printf("Error, not inMerse!n");return 0;if(isk != k)for(j=0; j<=n-1; j+)u=k*n+j; M=isk*n+j;p=au; au=aM; aM=

35、p; if(jsk!=k)for(i=0; i<=n-1; i+)u=i*n+k;M=i*n+jsk;p=au; au=aM; aM=p;l=k*n+k;al=1.0/al;for(j=0; j<=n-1; j+)if(j!=k)u=k*n+j; au=au*al;for(i=0; i<=n-1; i+)if(i!=k)for(j=0; j<=n-1; j+)if(j!=k)u=i*n+j;au=au-ai*n+k*ak*n+j; for(i=0; i<=n-1; i+)if(i!=k)u=i*n+k; au=-au*al;for(k=n-1; k>=0;

36、k-)if (jsk!=k)for(j=0; j<=n-1; j+) u=k*n+j; M=jsk*n+j;p=au; au=aM; aM=p; if(isk!=k)for(i=0; i<=n-1; i+) u=i*n+k; M=i*n+isk;p=au; au=aM; aM=p;delete is;delete js;return 1;/输出m行、n列的矩阵avoid shuchu(double *a, int m, int n)int i,j;for(i=0;i<m;i+)for(j=0;j<n;j+)printf("%9lf ",ai*n+j)

37、;printf("n");void main()zuobian();youbian();double a11,a12,a13,a21,a22,a23,b11,b12,b13,b21,b22,b23,c11,c12,c13,c21,c22,c23;double x1,y1,x2,y2,X1,Y1,Z1,X2,Y2,Z2,N1,N2,Xt,Yt,Zt,Bx,By,Bz;double f=0.15;scanf("x1=%lf,y1=%lf,x2=%lf,y2=%lf",&x1,&y1,&x2,&y2); /scanf("

38、;y1=%lf",&y1);/printf("x2="); /scanf("x2=%lf",&x2);/旋转矩阵a11=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01); a12=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);a13=(-1)*sin(p01)*cos(w01);a21=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02); a22=(-1)*cos(p02)*sin(k02)-sin

39、(p02)*sin(w02)*cos(k02);a23=(-1)*sin(p02)*cos(w02);b11=cos(w01)*sin(k01);b12=cos(w01)*cos(k01);b13=(-1)*sin(w01);b21=cos(w02)*sin(k02);b22=cos(w02)*cos(k02);b23=(-1)*sin(w02);c11=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);c12=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);c13=cos(p01)*cos(w01);c21=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);c22=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);c23=cos(p02)*cos(w02);/计算各待定点的像空间辅助坐标X1=a11*x1/100