中南大学本科生课程设计(实践)任务书、设计报告 (摄影测量与遥感概论)题目空间后方-前方交会学生姓名指导教师邹峥嵘学院地球科学与信息物理学院专业班级测绘0902班学生学号一、实验目的通过对数字影像空间后交前交的程序设计实验，要求我们进一步理解和掌握影像外方位元素的有关理论、原理和方法。

利用计算机程序设计语言编写摄影测量空间交会软件进行快速确定影像的外方位元素及其精度，然后通过求得的外方位元素求解未知点的地面摄影测量坐标，达到通过摄影测量量测地面地理数据的目的。

二、实验要求➢用C、VB或者Matlab编写空间后方交会-前方交会计算机程序。

➢提交实验报告：程序框图，程序源代码、计算结果及体会。

➢计算结果：地面点坐标、外方位元素及精度。

➢完成时间：2011年11月17日。

三、实验数据214.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，x0=0，y0=0四、实验思路➢利用后方交会得出两张像片各自的外方位元素1)获取已知数据：从摄影资料中插曲像片比例尺、平均航高、内方位元素以及控制点的地面摄影测量坐标及对应的像点坐标。

2)确定未知数的初始值：在竖直摄影的情况下，胶原素的初始值为0，线元素其中Zs=m*f+，Xs=，Ys=。

3)计算旋转矩阵R。

4)逐点计算像点坐标的近似值：利用共线方程。

5)组成误差方程并法化。

6)解求外方位元素。

7)检查计算是否收敛。

➢利用解求出的外方位元素进行前方交会1)用各自像片的角元素计算出左右像片的旋转矩阵R1和R2。

2)根据左右像片的外方位元素计算摄影基线分量Bx,By，Bz。

3)逐点计算像点的空间辅助坐标。

4) 计算投影系数。

5) 计算未知点的地面摄影测量坐标。

6) 重复以上步骤完成所有点的地面坐标的计算。

五、 实验过程➢ 程序流程框图不满足限差则重复计算➢程序中的主要函数设计子函数（矩阵求积multiply，计算函数Resection，矩阵转置transpose，矩阵求逆inMerse1，输出函数shuchu，左片的外方位元素求解函数zuobian。

右片的外方位元素求解函数youbian。

）➢程序源代码#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);//求矩阵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=fopen("F:\ground.txt","r")) == NULL){printf("Open file error!");return;}//像点坐标和地面点坐标double imagecontrol[4][2]={0.0};double groundcontrol[4][3]={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", &imagecontrol[i][j]);imagecontrol[i][j] /= 1000.0;}for(k=0; k<3; k++){fscanf(fp1, "%lf", &groundcontrol[i][k]);}}fclose(fp);fclose(fp1);//计算外方位元素初始值for( i=0;i<4;i++){Xs1+=groundcontrol[i][0];Ys1+=groundcontrol[i][1];Zs1+=groundcontrol[i][2];}Xs1/=4.0;Ys1/=4.0;Zs1/=4.0;Zs1+=m*f;double R[3][3]={0.0};double L3=0.0,L1=0.0,L2=0.0;double L[8][1]={0.0},x=0.0,y=0.0;double A[8][6]={0.0},A T[6][8]={0.0},ATA[6][6]={0.0},B[6][8]={0.0};double V[6][1]={0.0};int n=0;do{//计算旋转矩阵R[0][0]=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);R[0][1]=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R[0][2]=(-1)*sin(p01)*cos(w01);R[1][0]=cos(w01)*sin(k01);R[1][1]=cos(w01)*cos(k01);R[1][2]=(-1)*sin(w01);R[2][0]=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R[2][1]=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R[2][2]=cos(p01)*cos(w01);for(i=0,j=0;j<4;i+=2,j++){//计算像点坐标的近似值L1=R[0][0]*(groundcontrol[j][0]-Xs1)+R[1][0]*(groundcontrol[j][1]-Ys1)+R[2][0]*(ground control[j][2]-Zs1);L2=R[0][1]*(groundcontrol[j][0]-Xs1)+R[1][1]*(groundcontrol[j][1]-Ys1)+R[2][1]*(ground control[j][2]-Zs1);L3=R[0][2]*(groundcontrol[j][0]-Xs1)+R[1][2]*(groundcontrol[j][1]-Ys1)+R[2][2]*(ground control[j][2]-Zs1);x=(-1)*f*L1/L3;y=(-1)*f*L2/L3;//计算常数项L[2*j][0]=imagecontrol[j][0]-x;L[2*j+1][0]=imagecontrol[j][1]-y;//计算系数矩阵A[i][0]=(R[0][0]*f+R[0][2]*imagecontrol[j][0])/L3;A[i][1]=(R[1][0]*f+R[1][2]*imagecontrol[j][0])/L3;A[i][2]=(R[2][0]*f+R[2][2]*imagecontrol[j][0])/L3;A[i][3]=imagecontrol[j][1]*sin(w01)-((imagecontrol[j][0]/f)*(imagecontrol[j][0]*cos(k01)-i magecontrol[j][1]*sin(k01))+f*cos(k01))*cos(w01);A[i][4]=(-1)*f*sin(k01)-(imagecontrol[j][0]/f)*(imagecontrol[j][0]*sin(k01)+imagecontrol[j] [1]*cos(k01));A[i][5]=imagecontrol[j][1];A[i+1][0]=(R[0][1]*f+R[0][2]*imagecontrol[j][1])/L3;A[i+1][1]=(R[1][1]*f+R[1][2]*imagecontrol[j][1])/L3;A[i+1][2]=(R[2][1]*f+R[2][2]*imagecontrol[j][1])/L3;A[i+1][3]=(-1)*imagecontrol[j][0]*sin(w01)-((imagecontrol[j][1]/f)*(imagecontrol[j][0]*cos (k01)-imagecontrol[j][1]*sin(k01))-f*sin(k01))*cos(w01);A[i+1][4]=(-1)*f*cos(k01)-(imagecontrol[j][1]/f)*(imagecontrol[j][0]*sin(k01)+imagecontro l[j][1]*cos(k01));A[i+1][5]=(-1)*imagecontrol[j][0];}transpose(&A[0][0],&AT[0][0],8,6);multiply(&AT[0][0],&A[0][0],&ATA[0][0],6,8,6);inMerse1(*ATA,6);multiply(&ATA[0][0],&AT[0][0],&B[0][0],6,6,8);multiply(&B[0][0],&L[0][0],&V[0][0],6,8,1);Xs1+=V[0][0];Ys1+=V[1][0];Zs1+=V[2][0];p01+=V[3][0];w01+=V[4][0];k01+=V[5][0];n++;}while(fabs(V[3][0])>=0.00001||fabs(V[4][0])>=0.00001||fabs(V[5][0])>=0.00001);R[0][0]=cos(p01)*cos(k01)-sin(p01)*sin(w01)*sin(k01);R[0][1]=(-1)*cos(p01)*sin(k01)-sin(p01)*sin(w01)*cos(k01);R[0][2]=(-1)*sin(p01)*cos(w01);R[1][0]=cos(w01)*sin(k01);R[1][1]=cos(w01)*cos(k01);R[1][2]=(-1)*sin(w01);R[2][0]=sin(p01)*cos(k01)+cos(p01)*sin(w01)*sin(k01);R[2][1]=(-1)*sin(p01)*sin(k01)+cos(p01)*sin(w01)*cos(k01);R[2][2]=cos(p01)*cos(w01);//进行未知数的精度评定double A V[8][1],X[8][1],XT[1][8],XTX[1][1],mo,D[6][6],mi[6];multiply(&A[0][0],&V[0][0],&A V[0][0],8,6,1);for(i=0;i<8;i++)X[i][0]=A V[i][0]-L[i][0];transpose(&X[0][0],&XT[0][0],8,1);multiply(&XT[0][0],&X[0][0],&XTX[0][0],1,8,1);mo=XTX[0][0]/2;for(i=0;i<6;i++)for(j=0;j<6;j++)D[i][j]=mo*ATA[i][j];for(i=0;i<6;i++){mi[i]=sqrt(D[i][i]);}printf("左片结果为：\n\n");printf("旋转矩阵R为：\n");shuchu(&R[0][0],3,3);printf("外方为元素为：\n");printf("Xs1=%lf\n",Xs1);printf("Ys1=%lf\n",Ys1);printf("Zs1=%lf\n",Zs1);printf("p01=%lf\n",p01);printf("w01=%lf\n",w01);printf("k01=%lf\n",k01);printf("各外方位元素精度为：\n");printf("m1=%lf\nm2=%lf\nm3=%lf\nm4=%lf\nm5=%lf\nm6=%lf\n",mi[0],mi[1],mi[2],mi[3],mi[ 4],mi[5]);printf("迭代次数为：%d\n\n\n\n",n);fclose(fp);}void youbian(){FILE *fp = NULL;FILE *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 imagecontrol[4][2]={0.0};double groundcontrol[4][3]={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", &imagecontrol[i][j]);imagecontrol[i][j] /= 1000.0;}for(k=0; k<3; k++){fscanf(fp1, "%lf", &groundcontrol[i][k]);}}fclose(fp);fclose(fp1);//计算外方位元素初始值for( i=0;i<4;i++){Xs2+=groundcontrol[i][0];Ys2+=groundcontrol[i][1];Zs2+=groundcontrol[i][2];}Xs2/=4.0;Ys2/=4.0;Zs2/=4.0;Zs2+=m*f;double R[3][3]={0.0};double L3=0.0,L1=0.0,L2=0.0;double L[8][1]={0.0},x=0.0,y=0.0;double A[8][6]={0.0},A T[6][8]={0.0},ATA[6][6]={0.0},B[6][8]={0.0};double V[6][1]={0.0};int n=0;do{R[0][0]=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R[0][1]=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R[0][2]=(-1)*sin(p02)*cos(w02);R[1][0]=cos(w02)*sin(k02);R[1][1]=cos(w02)*cos(k02);R[1][2]=(-1)*sin(w02);R[2][0]=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R[2][1]=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R[2][2]=cos(p02)*cos(w02);for(i=0,j=0;j<4;i+=2,j++){//计算像点坐标的近似值L1=R[0][0]*(groundcontrol[j][0]-Xs2)+R[1][0]*(groundcontrol[j][1]-Ys2)+R[2][0]*(ground control[j][2]-Zs2);L2=R[0][1]*(groundcontrol[j][0]-Xs2)+R[1][1]*(groundcontrol[j][1]-Ys2)+R[2][1]*(ground control[j][2]-Zs2);L3=R[0][2]*(groundcontrol[j][0]-Xs2)+R[1][2]*(groundcontrol[j][1]-Ys2)+R[2][2]*(ground control[j][2]-Zs2);x=(-1)*f*L1/L3;y=(-1)*f*L2/L3;//计算常数项矩阵L[2*j][0]=imagecontrol[j][0]-x;L[2*j+1][0]=imagecontrol[j][1]-y;//计算系数矩阵A[i][0]=(R[0][0]*f+R[0][2]*imagecontrol[j][0])/L3;A[i][1]=(R[1][0]*f+R[1][2]*imagecontrol[j][0])/L3;A[i][2]=(R[2][0]*f+R[2][2]*imagecontrol[j][0])/L3;A[i][3]=imagecontrol[j][1]*sin(w02)-((imagecontrol[j][0]/f)*(imagecontrol[j][0]*cos(k02)-i magecontrol[j][1]*sin(k02))+f*cos(k02))*cos(w02);A[i][4]=(-1)*f*sin(k02)-(imagecontrol[j][0]/f)*(imagecontrol[j][0]*sin(k02)+imagecontrol[j] [1]*cos(k02));A[i][5]=imagecontrol[j][1];A[i+1][0]=(R[0][1]*f+R[0][2]*imagecontrol[j][1])/L3;A[i+1][1]=(R[1][1]*f+R[1][2]*imagecontrol[j][1])/L3;A[i+1][2]=(R[2][1]*f+R[2][2]*imagecontrol[j][1])/L3;A[i+1][3]=(-1)*imagecontrol[j][0]*sin(w02)-((imagecontrol[j][1]/f)*(imagecontrol[j][0]*cos (k02)-imagecontrol[j][1]*sin(k02))-f*sin(k02))*cos(w02);A[i+1][4]=(-1)*f*cos(k02)-(imagecontrol[j][1]/f)*(imagecontrol[j][0]*sin(k02)+imagecontro l[j][1]*cos(k02));A[i+1][5]=(-1)*imagecontrol[j][0];}transpose(&A[0][0],&AT[0][0],8,6);multiply(&AT[0][0],&A[0][0],&ATA[0][0],6,8,6);inMerse1(*ATA,6);multiply(&ATA[0][0],&AT[0][0],&B[0][0],6,6,8);multiply(&B[0][0],&L[0][0],&V[0][0],6,8,1);Xs2+=V[0][0];Ys2+=V[1][0];Zs2+=V[2][0];p02+=V[3][0];w02+=V[4][0];k02+=V[5][0];n++;}while(fabs(V[3][0])>=0.00001||fabs(V[4][0])>=0.00001||fabs(V[5][0])>=0.00001);R[0][0]=cos(p02)*cos(k02)-sin(p02)*sin(w02)*sin(k02);R[0][1]=(-1)*cos(p02)*sin(k02)-sin(p02)*sin(w02)*cos(k02);R[0][2]=(-1)*sin(p02)*cos(w02);R[1][0]=cos(w02)*sin(k02);R[1][1]=cos(w02)*cos(k02);R[1][2]=(-1)*sin(w02);R[2][0]=sin(p02)*cos(k02)+cos(p02)*sin(w02)*sin(k02);R[2][1]=(-1)*sin(p02)*sin(k02)+cos(p02)*sin(w02)*cos(k02);R[2][2]=cos(p02)*cos(w02);//进行未知数的精度评定double A V[8][1],X[8][1],XT[1][8],XTX[1][1],mo,D[6][6],mi[6];multiply(&A[0][0],&V[0][0],&A V[0][0],8,6,1);for(i=0;i<8;i++)X[i][0]=A V[i][0]-L[i][0];transpose(&X[0][0],&XT[0][0],8,1);multiply(&XT[0][0],&X[0][0],&XTX[0][0],1,8,1);mo=XTX[0][0]/2;for(i=0;i<6;i++)for(j=0;j<6;j++)D[i][j]=mo*ATA[i][j];for(i=0;i<6;i++){mi[i]=sqrt(D[i][i]);}printf("右片结果为：\n\n");printf("旋转矩阵R为：\n");shuchu(&R[0][0],3,3);printf("外方为元素为：\n");printf("Xs2=%lf\n",Xs2);printf("Ys2=%lf\n",Ys2);printf("Zs2=%lf\n",Zs2);printf("p02=%lf\n",p02);printf("w02=%lf\n",w02);printf("k02=%lf\n",k02);printf("各外方位元素精度为：\n");printf("m1=%lf\nm2=%lf\nm3=%lf\nm4=%lf\nm5=%lf\nm6=%lf\n",mi[0],mi[1],mi[2],mi[3],mi[ 4],mi[5]);printf("迭代次数为：%d\n",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++){b[j*m+i] = a[i*n+j];}}}//矩阵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 += a[i*n+k]*b[k*l+j];}c[i*l+j]=t;}}}//求矩阵a的逆int inMerse1(double *a, int n){int *is, *js, i, j, k, l, u, M;double d,p;is = new int[n];js = new int[n];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(a[l]);if (p>d) { d=p; is[k]=i; js[k]=j;}}}if(d+1.0==1.0){delete []is;delete []js;printf("Error, not inMerse!\n");return 0;}if(is[k] != k)for(j=0; j<=n-1; j++){u=k*n+j;M=is[k]*n+j;p=a[u];a[u]=a[M];a[M]=p;}if(js[k]!=k)for(i=0; i<=n-1; i++){u=i*n+k;M=i*n+js[k];p=a[u];a[u]=a[M];a[M]=p;}l=k*n+k;a[l]=1.0/a[l];for(j=0; j<=n-1; j++)if(j!=k){u=k*n+j; a[u]=a[u]*a[l];}for(i=0; i<=n-1; i++)if(i!=k)for(j=0; j<=n-1; j++)if(j!=k){u=i*n+j;a[u]=a[u]-a[i*n+k]*a[k*n+j];}for(i=0; i<=n-1; i++)if(i!=k){u=i*n+k;a[u]=-a[u]*a[l];}}for(k=n-1; k>=0; k--){if (js[k]!=k){for(j=0; j<=n-1; j++){u=k*n+j; M=js[k]*n+j;p=a[u]; a[u]=a[M]; a[M]=p;}}if(is[k]!=k){for(i=0; i<=n-1; i++){u=i*n+k; M=i*n+is[k];p=a[u]; a[u]=a[M]; a[M]=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 ",a[i*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;double 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=(-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(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/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)/(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-1500;printf("Xt=%lf\nYt=%lf\nZt=%lf\n",Xt,Yt,Zt);}➢实验结果左片右片地面摄影测量坐标点号x y x y X Y Z GCP116.01279.963-73.9378.7065083.2055852.099527.925 GCP288.5681.134-5.25278.1845780.025906.365571.549 GCP313.362-79.37-79.122-78.8795210.8794258.446461.81 GCP482.24-80.027-9.887-80.0895909.2644314.283455.484 151.75880.555-39.95378.4635431.489 5879.367 549.723 214.618-0.231-76.0060.0365147.388 5055.549 485.001 349.88-0.782-42.201-1.0225495.786 5082.726 506.677 486.14-1.346-7.706-2.1125844.173 5109.861 528.420 548.035-79.962-44.438-79.7365559.940 4268.185 463.523六、实验总结从编写思路及步骤，然后是梳理框架，要用到的子函数，子函数参数，怎么编写，有不懂的东西网上参考资料；再是自己动手编程，发现编程要求很细腻，出不得一点差错，程序编写出来后，有很多的错误，要求逐一进行改正，修改；最后才得以运行。