作业4--空间后方交会资料.docx

摄影测量作业报告之空间后方交会陈闻亚20080729作业报告空间后方交会专 业：测绘工程班 级： 2008 级班姓 名：陈闻亚指导教师：陈强2010年 4月 16日1作业任务 32作业思想 33作业条件及数据 34作业过程 35源程序 46计算结果 177心得体会与建议 17-1 -摄影测量作业报告之空间后方交会陈闻亚200807291作业任务计算近似垂直摄影情况下后方交会解。即利用摄影测量空间后方交会的方法，获取相 片的6个外方位元素。限差为 0.1。2作业思想利用摄影测量空间后方交会的方法求解。该方法的基本思想是利用至少三个一直地面 控制点的坐标 A (Xa, Ya, Za)、B (Xb, Yb, Zb) C (Xc, Yc, Zc),与其影像上对应的 三个像点的影像坐标a (%,ya)、b(xb,yb)、c(xc,yc),根据共线方程，反求该相片的外方位元素 Xs、Ys、Zs、6 3、ko3作业条件及数据已知摄影机主距f=153.24mm ,四对点的像点坐标与相应的地面坐标列入下表:表1点万像点坐标地面坐标x (mm)y(mm)X(m)Y(m)Z(m)1-86.15-68.9936589.4125273.322195.172-53.40-82.2137631.0831324.51.728.69314.7876.6339100.9724934.982386.50410.4664.4340426.5430319.81757.314作业过程4. 1获取已知数据相片比例尺1/m=1 : 10000,内方位元素f=153.24mm , x, y;获取控制点的地面测量 坐标 Xt、Yt、Zt。4. 2量测控制点的像点坐标：本次作业中为已知。见表 1。4. 3确定未知数的初始值：在近似垂直摄影情况下，月射M素的初始值为0,即 而=助=咫=0；线元素中，Zs0=H=mf=1532.4m , Xs。、Ys0的取值可用四个控制点坐标的平均值，即：Xs0= 144、Xi =38437.00-9 -1丫 so=44、丫 =89106.62i 14. 4计算旋转矩阵R:利用胶原素的近似值计算方向余弦值，组成R阵。4. 5逐点计算像点坐标的近似值：x) (y)。利用未知数的近似值按共线方程式计算控制点像点坐标的近似值(4. 6组成误差方程:陈闻亚20080729逐点计算误差方程式的系数和常数项。4. 7组成法方程式：计算法方程的系数矩阵ATA与常数项ATL。4. 8求解外方位元素：根据法方程，由式 X= (AtA) -1 ATL解求外方位元素改正数，并与相应的近似值求和， 得到外方位元素新的近似值。4. 9求解外方位元素：将求得的外方位元素的改正数与规定的限差(0.1)比较，小于限差则计算终止，否则用新的近似值重复第 4.4至4.8步骤的计算，知道满足要求为止。5源程序#include <stdio.h>#include <stdlib.h>#include <math.h> const double PRECISION=1e-5;typedef double DOUBLE5;int InputData(int &Num, DOUBLE *&Data,double &m,double &f);int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f); int InverseMatrix(double *matrix,const int &row);int main(int argc, char* argv) DOUBLE *Data=NULL;int Num;double f(0),m(0);if(InputData(Num,Data,m,f)if (Data!=NULL) delete 口Data; return 1;if(Resection(Num,Data,m,f)if (Data!=NULL) delete 口Data;return 1;陈闻亚20080729if (Data!=NULL)delete 口Data;printf("解算完毕n");doprintf("计算结果保存于 "结果.txt"文件中n"”请选择操作(输入P打开结果数据，R打开原始数据，其它退出程序)：”);fflush(stdin);刷新输入流char order=getchar();if (P=order | p=order)system("结果.txt");else if (R=order | r=order)system("data.txt");elsebreak;system("cls");while(1);system("PAUSE");return 0;/* 函数名：InputData* 函数介绍：从文件(data.txt)中读取数据，* 文件格式如下：* 点数m(未知写作0)* 内方位元素(f x0 y0)* 编号x y X Y Z* 实例:4 0153.24 0 01 -86.15-68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100.97 24934.98 2386.504 10.46 64.43 40426.54 30319.81 757.31 *参数:(in/out)Num(点数),陈闻亚20080729*(in/out)Data(存放数据),m,f,x0,y0*返回值：int , 0成功，1文件打开失败，2控制点个*数不足，3文件格式错误*/int InputData(int &Num, DOUBLE *&Data,double &m,double &f) double x0,y0;FILE *fp_input;if (!(fp_input=fopen("data.txt","r")return 1;fscanf(fp_input,"%d%lf",&Num,&m);if (Num<4)return 2;fscanf(fp_input,"%lf%lf%lf,&f,&x0,&y0);f/=1000;if (m<0 | f<0)return 3;Data=new DOUBLENum;double *temp= new doubleNum-1;double scale=0;int i;for (i=0;i<Num;i+)读取数据，忽略编号if(fscanf(fp_input,"%*d%lf%lf%lf%lf%lf", &Datai0,&Datai1,&Datai2, &Datai3,&Datai4)!=5)return 3;单位换算成mDatai0/=1000.0;Datai1/=1000.0;陈闻亚20080729/如果m未知则归算其值if (0=m)for (i=0;i<Num-1;i+)tempi=(Datai2-Datai+12)/(Datai0-Datai+10)+ (Datai3-Datai+13)/(Datai1-Datai+11);scale+=tempi/2.0;m=scale/(Num-1);fclose(fp_input);delete 口temp;return 0;/* 函数名：MatrixMul* 函数介绍：求两个矩阵的积，* 参数：Jz1(第一个矩阵工row(第一个矩阵行数)，* Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个* 矩阵列数),(out)JgJz(存放结果矩阵)*返回值：void*/void MatrixMul(double *Jz1,const int &row,double *Jz2, const int &line,const int &com,double *JgJz)for (int i=0;i<row;i+)for (int j=0;j<line;j+)double temp=0;for (int k=0;k<com;k+)temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j);*(JgJz+i*line+j尸temp;*陈闻亚20080729* 函数名：OutPut* 函数介绍：向结果.txt文件输出数据* 参数：Q协因数阵，m精度，m0单位权中误差，6个外* 方位元素，旋转矩阵*返回值：int, 0成功，1失败*/int OutPut(const double *&Q,const double *&m,const double &m0, const double &Xs,const double &Ys,const double &Zs, const double &Phi,const double &Omega, const double &Kappa,const double *R)FILE *fp_out;if (!(fp_out=fopen("结果.txt","w")return 1;FILE *fp_input;if (!(fp_input=fopen("data.txt","r")return 1; fprintf(fp_out, "*""*""*""*n");fprintf(fp_out,"n空间后方交会程序(CC+)n测绘一班n" ”学号：20080729n姓名：陈闻亚nn");fprintf(fp_out,"*""*""*""*n");fprintf(fp_out,"已知数据：nn 已知点数：");int num;double temp,x,y;fscanf(fp_input,"%d%lf",&num,&temp);fprintf(fp_out,"%dn",num);fprintf(fp_out,"摄影比例尺(0表示其值位置)：");fprintf(fp_out,"%10.0lfn",temp);fprintf(fp_out,"内方位元素(f x0 y0):"); fscanf(fp_input,"%lf%lf%lf",&temp,&x,&y);陈闻亚20080729fprintf(fp_out,"%10lft%10lft%10lfn",temp,x,y);for (int i=0;i<num;i+)double temp5;fscanf(fp_input,"%*d%lf%lf%lf%lf%lf", &temp0,&temp1,&temp2,&temp3,&temp4);fprintf(fp_out,"%3dt%10lft%10lft%10lft%10lft%10lfn", i+1,temp0,temp1,temp2,temp3,temp4);fclose(fp_input);fprintf(fp_out,"*""*""*""*n，);fprintf(fp_out,"计算结果如下：nn外方位元素：n");fprintf(fp_out,"tXs=%10lfn",Xs);fprintf(fp_out,"tYs=%10lfn",Ys);fprintf(fp_out,"tZs=%10lfn",Zs);fprintf(fp_out,"tPhi=%10lfn",Phi);fprintf(fp_out,"tOmega=%10lfn",Omega);fprintf(fp_out,"tKappa=%101fnn",Kappa);fprintf(fp_out,"旋转矩阵：n");for (i=0;i<3;i+)fprintf(fp_out,"t");for (int j=0;j<3;j+)fprintf(fp_out,"%10lft",*(R+i*3+j);fprintf(fp_out,"n");fprintf(fp_out,"n 单位权中误差：%10lfnn",m0);fprintf(fp_out,"协因数阵：n");for (i=0;i<6;i+)fprintf(fp_out,"t");for (int j=0;j<6;j+) 陈闻亚20080729fprintf(fp_out,"%20lft",*(Q+i*6+j);fprintf(fp_out,"n");fprintf(fp_out,"n 外方位元素精度：")；for (i=0;i<6;i+)fprintf(fp_out,"%10lft",mi);fprintf(fp_out,"n");fprintf(fp_out, "*""*"”*”"*n");fclose(fp_out);return 0; /*函数名：Resection函数介绍：计算参数：Num（点数）,Data（数据），m,f（焦距）,x0,y0返回值：int, 0成功，其它失败 */int Resection(const int &Num,const DOUBLE *&Data,const double &m, const double &f)double Xs=0,Ys=0,Zs=0; int i,j;/设置初始值for (i=0;i<Num;i+)Xs+=Datai2;Ys+=Datai3;Xs/=Num;Ys/=Num;Zs=m*f;double Phi(0),Omega(0),Kappa(0);摄影测量作业报告之空间后方交会陈闻亚20080729double R33=0.0;double *L=new double2*Num;typedef double Double66;Double6 *A=new Double62*Num;double *AT=new double2*Num*6;double *ATA=new double6*6;double *ATL=new double6;double *Xg=new double6;/迭代计算do旋转矩阵R00=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa);R01=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa);R02=-sin(Phi)*cos(Omega);R10=cos(Omega)*sin(Kappa);R11=cos(Omega)*cos(Kappa);R12=-sin(Omega);R20=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa);R21=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa);R22=cos(Phi)*cos(Omega);for (i=0;i<Num;i+)double X=R00*(Datai2-Xs)+R10*(Datai3-Ys)+ R20*(Datai4-Zs);double Y=R01*(Datai2卜Xs)+R11*(Datai3-Ys)+ R21*(Datai4-Zs);double Z=R02*(Datai2-Xs)+R12*(Datai3-Ys)+ R22*(Datai4-Zs);double xxx,yyy;xxx=-f*X/Z;yyy=-f*Y/Z;常数项L2*i=Datai0-(-f*X/Z);L2*i+1=Datai1-(-f*Y/Z);A2*i0=(R00*f+R02*(xxx)/Z;A2*i1=(R10*f+R12*(xxx)/Z;A2*i2=(R20*f+R22*(xxx)/Z;-11 -摄影测量作业报告之空间后方交会陈闻亚20080729A2*i3=(yyy)*sin(Omega)-(xxx)/f)* (xxx)*cos(Kappa)-(yyy)*sin(Kappa)+ f*cos(Kappa)*cos(Omega);A2*i4=-f*sin(Kappa)-(xxx)/f)*(xxx)* sin(Kappa)+(yyy)*cos(Kappa);A2*i5=(yyy);A2*i+10=(R01*f+R02*(yyy)/Z;A2*i+11=(R11*f+R12*(yyy)/Z;A2*i+12=(R21*f+R22*(yyy)/Z;A2*i+13=-(xxx)*sin(Omega)-(yyy)/f)* (xxx)*cos(Kappa)-(yyy)*sin(Kappa)- f*sin(Kappa)*cos(Omega);A2*i+14=-f*cos(Kappa)-(yyy)/f)*(xxx)* sin(Kappa)+(yyy)*cos(Kappa);A2*i+15=-(xxx);求矩阵A的转置矩阵ATfor (i=0;i<2*Num;i+)for (j=0;j<6;j+)*(AT+j*2*Num+i)=Aij;求 ATAMatrixMul(AT,6,&A00,6,2*Num,ATA);if(InverseMatrix(ATA,6) return 1;MatrixMul(AT,6,L,1,2*Num,ATL);MatrixMul(ATA,6,A TL,1,6,Xg);Xs+=Xg0;Ys+=Xg1;Zs+=Xg2;Phi+=Xg3;Omega+=Xg4;Kappa+=Xg5; while(fabs(Xg0)>=PRECISION |fabs(Xg1)>=PRECISION |fabs(Xg2)>=PRECISION |fabs(Xg3)>=PRECISION | fabs(Xg4)>=PRECISION | (Xg5)>=PRECISION);陈闻亚20080729/注：协因数阵，旋转矩阵等计算本应该使用最后外方位元素值,由于变换很小忽略double *Q=ATA;double *V=new double2*Num;MatrixMul(&A00,2*Num,Xg,1,6,V);double VTV=0;for(i=0;i<2*Num;i+)Vi-=Li;VTV+=Vi*Vi;double m0=sqrt(VTV/(2*Num-6);double *mm=new double6;for (i=0;i<6;i+)mmi=sqrt(*(Q+i*6+i)*m0;OutPut(Q,mm,m0,Xs,Ys,Zs,Phi,Omega,Kappa,&R00);delete 口L;delete 口A;delete 口AT;delete 口ATA;delete 口ATL;delete Xg;delete 口mm;delete 口V;return 0;void swap(double &a,double &b)double temp=a;a=b;b=temp;/* 函数名：InverseMatrix*函数介绍：求矩阵的逆(高斯-约当法)*输入参数：（in/out）matrix （矩阵首地址）陈闻亚20080729* (in) row (矩阵阶数)* 输出参数：matrix (原矩阵的逆矩阵)* 返回值：int , 0成功，1失败* 调用函数：swap(double&,double&)*/int InverseMatrix(double *matrix,const int &row)double *m=new doublerow*row;double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;i<row;i+)for (j=0;j<row;j+)*pt=*ptemp;ptemp+;pt+;int k;int *is=new introw,*js=new introw;for (k=0;k<row;k+)double max=0;/全选主元寻找最大元素for (i=k;i<row;i+)for (j=k;j<row;j+)if (fabs(*(m+i*row+j)>max)max=*(m+i*row+j);isk=i;jsk=j;-13 -摄影测量作业报告之空间后方交会陈闻亚20080729if (0 = max)(return 1;)行交换if (isk!=k)(for (i=0;i<row;i+)(swap(*(m+k*row+i),*(m+isk*row+i);)列交换if (jsk!=k)(for (i=0;i<row;i+)(swap(*(m+i*row+k),*(m+i*row+jsk);)*(m+k*row+k)=1/(*(m+k*row+k);for (j=0;j<row;j+)(if (j!=k)(*(m+k*row+j)*=*(m+k*row+k);)for (i=0;i<row;i+)(if (i!=k)(for (j=0;j<row;j+)(if(j!=k)(*(m+i*row+j)-=*(m+i*row+k)*(m+k*row+j);陈闻亚20080729)for (i=0;i<row;i+)(if(i!=k)(*(m+i*row+k)*=-(*(m+k*row+k);)int r;恢复行列for (r=row-1;r>=0;r-)(if (jsr!=r)(for (j=0;j<row;j+)(swap(*(m+r*row+j),*(m+jsr*row+j);)if (isr!=r)(for (i=0;i<row;i+)(swap(*(m+i*row+r),*(m+i*row+isr);)ptemp=matrix;pt=m;for (i=0;i<row;i+)(for (j=0;j<row;j+)(*ptemp=*pt;ptemp+;pt+;)delete 口is; delete js; delete 口m;return 0;6计算结果外方位元素：Xs=39795.452297旋转矩阵：0.997709-0.067526-0.004121Ys=27476.462210 Zs=7572.685927 （）=-0.003987 0.002114 行-0.0675780.0675340.0039870.997715-0.0021140.0018400.999990单位权中误差：0.0000077心得体会与建议有挑战才会有激情，成功了便更有满足感。此次作业，便是又一次的挑战。总体来说， 空间后方交会的知识，在学习时就有点懵，后来自己下来又看书加深了理解，把它弄懂了。 而编程方面所面临的问题就比较多了，上一学期就有过一个类似的作业，当时没有攻破编程这一题，便用的Excel解决，但事后也没有继续研究，实为遗憾。直到遇到这次作业，又是 同样的问题，才后悔当初没有继续学习把它解决了。其实主要总结起来就一个比较难的点： 编程实现矩阵的运算尤其是矩阵求逆。搜集资料时，在网上搜索到了一些矩阵求逆的代码， 但是有些也不能实现，找到合适的后，也要考虑把它运用到主程序的问题。不断调试，不断改进，这个过程还是有些难度的，最后还实现了用文本导入数据计算及文本输出结果。此次作业虽然比较波折， 但是收获还是很大， 也解决了一个可以算是拖延了很久的问题。虽然有些地方是效仿的，不过也是自己要弄懂以后，再试着自己一步步理出来。 对自己可以说是又一次提高。-17 -