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

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

2、序源代码、计算结果及体会。? 计算结果：地面点坐标、外方位元素及精度。? 完成时间： 2011 年 11 月 17 日。三、实验数据左片右片地面摄影测量坐标点号xyxyXYZGCP116.01279.963-73.9378.7065083.2055852.099527.925GCP288.5681.134-5.25278.1845780.025906.365571.549GCP313.362-79.37-79.122-78.8795210.8794258.446461.81GCP482.24-80.027-9.887-80.0895909.2644314.283455.484151.75880

3、.555-39.95378.463214.618-0.231-76.0060.036349.88-0.782-42.201-1.022486.14-1.346-7.706-2.112548.035-79.962-44.438-79.736f=150.000mm，x 0 =0 ，y 0=0四、实验思路? 利用后方交会得出两张像片各自的外方位元素1) 获取已知数据：从摄影资料中插曲像片比例尺、平均航高、内方位元素以及控制点的地面摄影测量坐标及对应的像点坐标。2) 确定未知数的初始值： 在竖直摄影的情况下， 胶原素的初始值为 0，线元素其中 Zs=m*f+，Xs=，Ys=。3) 计算旋转矩阵 R。4

4、) 逐点计算像点坐标的近似值：利用共线方程。5) 组成误差方程并法化。6) 解求外方位元素。7) 检查计算是否收敛。? 利用解求出的外方位元素进行前方交会1) 用各自像片的角元素计算出左右像片的旋转矩阵R1 和 R2。2) 根据左右像片的外方位元素计算摄影基线分量 Bx,By ，Bz。3) 逐点计算像点的空间辅助坐标。4) 计算投影系数。5) 计算未知点的地面摄影测量坐标。6) 重复以上步骤完成所有点的地面坐标的计算。五、实验过程? 程序流程框图后方交会函数确定 已知数据比例尺 m确定各 外方位元素初始值计算旋转矩阵逐点计算像点坐标不近似值满足限逐点计算误差方程系差数项，组成误差系数则矩阵 A

5、重复利用矩阵运算求解外计方位元素算检查是否满足限差若满足则输出外方位元素将整个过程作为一个函数继续进行右片的外方位元素求解求解各外方位元素精度空间前方交会利用已求得的角元素计算2 张像片各自的旋转矩阵利用已求得的线元素Xs1， Ys1， Zs1，p0，,w01，k01；Xs2 ，Ys2，Zs2，p02，w02， k02，计算基线分量：Bx=Xs2-Xs1 ；By=Ys2-Ys1 ；Bz=Zs2-Zs1 ；输入像片坐标，利用旋转矩阵求解想空间辅助坐标计算点投影系数：N1=(Bx*Z2-Bz*X2)/(X1*Z2-X2*Z1)；N2=(Bx*Z1-Bz*X1)/(X1*Z2-X2*Z1)；计算地面

6、摄影测量坐标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 。）? 程序源代码#include stdio.h#include math.hdouble Xs1,Xs2,Ys1,Ys2,Zs1,Zs2,p

7、01,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);/求矩阵 a 的逆int inMerse1(double *a, int n);/输出 m 行、 n 列的矩阵avoid shuchu(double *a, int m

8、, 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=fopen(F:ground.txt,r) = NULL)printf(Open file error!);return;/像点坐标和地面点坐标double imagecontrol42=0.0;double groundco

9、ntrol43=0.0;/摄影比例尺分母double m = 9943;double f=0.15;long i,j,k;for(i=0; i4; i+)for(j=0; j2; j+)fscanf(fp, %lf, &imagecontrolij);imagecontrolij /= 1000.0;for(k=0; k3; k+)fscanf(fp1, %lf, &groundcontrolik);fclose(fp);fclose(fp1);/计算外方位元素初始值for( i=0;i4;i+)Xs1+=groundcontroli0;Ys1+=groundcontroli1;Zs1+=gr

10、oundcontroli2;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 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;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

11、=(-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);for(i=0,j=0;j=0.00001|fabs(V40)=0.00001|fabs(V50)=0.00001);R00=cos(p01)*cos(k01)-sin(p01)

12、*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);/进行未知数的精度评定do

13、uble AV81,X81,XT18,XTX11,mo,D66,mi6;multiply(&A00,&V00,&AV00,8,6,1);for(i=0;i8;i+)Xi0=A Vi0-Li0;transpose(&X00,&XT00,8,1);multiply(&XT00,&X00,&XTX00,1,8,1);mo=XTX00/2;for(i=0;i6;i+)for(j=0;j6;j+)Dij=mo*ATAij;for(i=0;i6;i+)mii=sqrt(Dii);printf( 左片结果为：nn);printf( 旋转矩阵R 为： n);shuchu(&R00,3,3);printf( 外

14、方为元素为：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( 各外方位元素精度为：n);printf(m1=%lfnm2=%lfnm3=%lfnm4=%lfnm5=%lfnm6=%lfn,mi0,mi1,mi2,mi3,mi 4,mi5);printf( 迭代次数为：%dnnnn,n);fclose(fp);void youbian()FILE *fp = NULL;FIL

15、E *fp1 = NULL;if(fp=fopen(F:image2.txt,r) = NULL)printf(Open file error!);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; i4; i+)for(j=0; j2; j+)fscan

16、f(fp, %lf, &imagecontrolij);imagecontrolij /= 1000.0;for(k=0; k3; k+)fscanf(fp1, %lf, &groundcontrolik);fclose(fp);fclose(fp1);/计算外方位元素初始值for( i=0;i4;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=0.0,L2=0.0;doubl

17、e 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=sin(p02)*cos(k0

18、2)+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=0.00001|fabs(V40)=0.00001|fabs(V50)=0.00001);R00=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=c

19、os(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(w02);/进行未知数的精度评定double AV81,X81,XT18,XTX11,mo,D66,mi6;multiply(&A00,&V00,&AV00,8,6,1);for(i=0;i8;i+)Xi0=A Vi0-Li0;transpose

20、(&X00,&XT00,8,1);multiply(&XT00,&X00,&XTX00,1,8,1);mo=XTX00/2;for(i=0;i6;i+)for(j=0;j6;j+)Dij=mo*ATAij;for(i=0;i6;i+)mii=sqrt(Dii);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=%

21、lfn,w02);printf(k02=%lfn,k02);printf( 各外方位元素精度为：n);printf(m1=%lfnm2=%lfnm3=%lfnm4=%lfnm5=%lfnm6=%lfn,mi0,mi1,mi2,mi3,mi 4,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;im;i+)for(j=0;jn;j+)bj*m+i = ai*n+j;/矩阵 a 乘以矩阵 b

22、，结果存储在 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;im;i+)for(j=0;jl;j+)t=0;for(k=0;kn;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 = new intn;js = new

23、intn;for(k=0; k=n-1; k+)d=0.0;for(i=k; i=n-1; i+)for(j=k; jd) 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=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

24、(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=0; 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

25、 1;/输出 m 行、 n 列的矩阵avoid shuchu(double *a, int m, int n)int i,j;for(i=0;im;i+)for(j=0;jn;j+)printf(%9lf ,ai*n+j);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; doubl

26、e f=0.15;scanf(x1=%lf,y1=%lf,x2=%lf,y2=%lf,&x1,&y1,&x2,&y2);/scanf(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=(

27、-1)*cos(p02)*sin(k02)-sin(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(k0

28、1);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/1000+a12*y1/1000-a13*f;Y1=b11*x1/1000+b12*y1/1000-b13*f;Z1=c11*x1/1000+c12*y1/1000-c13*f;X2=a21*x2/1000+a22*y2/1000-a23*f;Y2=b21*x2/1000+b22*y2/1000-b23*f;Z2=c21*x2/1000+c22*y2/1000-c23*f;/计算摄影基线分量Bx=Xs2-Xs1;By=Ys2-Ys1;Bz=Zs2-Zs1;/计算点投影系数N1=(Bx*Z2-Bz*X2)/(X1*Z2-X2*Z1);N2=(Bx*Z1-Bz*X1