相对定向--课程设计报告

课程设计报告相对定向元素计算2010年5月15日1设计任务------------------------------------------------------------------------------------32作业原理---------------------------------------------------------------------------------------33已知条件及数据--------------------------------------------------------------------34作业过程---------------------------------------------------------------------------35程序流程-----------------------------------------------------------------------------46源程序-----------------------------------------------------------------------------57计算结果---------------------------------------------------------------------------108心得体会与建议-----------------------------------------------------------------------------101设计任务采用连续像对相对定向方法，计算出两张像片的相对定向元素，并运用一种高级语言解决此问题。2作业原理解析法相对定向原理：根据同名光线对对相交这一立体像对内在的几何关系，通过量测的像点坐标，用解析计算的方法解求相对定向元素，建立与地面相似的立体模型，确定模型点的三维坐标。相对定向的共面条件：B·（S1a1·S2a2）=0，即F=322111ZYXZYXBBBxZY连续像对的相对定向：连续像对法相对定向是以左像片为基准，求出右像片相对于左像片的五个定向元素222,,,,ZYBB.为了统一单位，吧bY，bZ两个基线元素改为角度形式表示，如下，μ和ν为极限的偏角和倾角。·btan·cosb·btan·XZXXXYbbb将上式代入共面条件方程，B=m·b，F=222111ZYXZYXbbbXXX=Xb2221111ZYXZYX=03已知条件及数据在一个航带内相邻两张像片上分别量测了6对同名点的像平面坐标，f=24mm，数据如下：点号左片右片x(mm)y(mm)x(mm)y(mm)1-6.0911.983-5.564-3.20227.0980.9247.694-2.83034.5381.0685.098-2.87846.8581.2087.429-2.5785-10.050-0.514-9.152-5.6426-8.0891.293-7.441-3.9814作业过程4．1原始数据的输入及像点坐标的预处理把在像点坐标量测仪器上量测出的各定向点的坐标换算成以像主点为原点的像点坐标并作系统误差改正。4．2确定相对定向元素的初始值：对于连续像对法进行相对定向，如果是航线的第一个像对，则左像片的角元素为零；若是后续像对，则由前一像对的相对定向算得本像对左像片的角元素。右片的三个角元素的初始值取零，基线分量by、bz也取零，bx取标准点位第一点的左右视差，即bx=(x1-x2)1。4．3计算左片的方向余弦值，组成旋转矩阵R1，计算左片各像点的像空间辅助坐标系（X1,Y1,Z1）。4．4计算右片的方向余弦值，组成旋转矩阵R2，计算基线分量by和bz。4．5计算右片各像点的空间辅助坐标（X2,Y2,Z2），计算各项点的投影系数N1，N2和上下视差Q；4．6逐点组成误差方程并法化，完成法方程系数矩阵和常数项矩阵的计算。4．7解法方程，求出相对定向元素的改正数。4．8计算相对定向元素的新值。根据法方程，由式X=（AtA）-1ATL解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。4．9检查所有的改正数是否小于限值0.00003rad，如满足条件，则结束相对定向计算。否则重复4.4~4.9。5程序流程6源程序#includeiostream.h#includefstream.h#includestdlib.h#includeiomanip.h#includemath.hconstintN=6;inti,j,k;doubleb,x[6][3],y[6][3],z[6][3]={0};//求转置矩阵templatetypenameT1,typenameT2voidTranspose(T1*mat1,T2*mat2,inta,intb){for(i=0;ia;i++)for(j=0;jb;j++)mat2[j][i]=mat1[i][j];return;}//求矩阵的乘积templatetypenameT1,typenameT2voidArray_mul(T1*mat1,T2*mat2,T2*result,inta,intb,intc){inti,j,k;for(i=0;ia;i++){for(j=0;jc;j++){result[i][j]=0;for(k=0;kb;k++)result[i][j]+=mat1[i][k]*mat2[k][j];}}return;}//求逆矩阵voidswap(double*a,double*b){doublec;c=*a;*a=*b;*b=c;};Inverse(doubleA[5][5],intn){inti,j,k;doubled;intJS[N],IS[N];for(k=0;kn;k++){d=0;for(i=k;in;i++)for(j=k;jn;j++){if(fabs(A[i][j])d){d=fabs(A[i][j]);IS[k]=i;JS[k]=j;};};if(d+1.0==1.0)return0;if(IS[k]!=k)for(j=0;jn;j++)swap(&A[k][j],&A[IS[k]][j]);if(JS[k]!=k)for(i=0;in;i++)swap(&A[i][k],&A[i][JS[k]]);A[k][k]=1/A[k][k];for(j=0;jn;j++)if(j!=k)A[k][j]=A[k][j]*A[k][k];for(i=0;in;i++)if(i!=k)for(j=0;jn;j++)if(j!=k)A[i][j]=A[i][j]-A[i][k]*A[k][j];for(i=0;in;i++)if(i!=k)A[i][k]=-A[i][k]*A[k][k];};for(k=n-1;k=0;k--){for(j=0;jn;j++)if(JS[k]!=k)swap(&A[k][j],&A[JS[k]][j]);for(i=0;in;i++)if(IS[k]!=k)swap(&A[i][k],&A[i][IS[k]]);};return1;}//原始数据导入voidInput(){doublem;j=0;i=0;ifstreamf1(左片各点像空间坐标.txt);if(!f1){cerr左片各点像空间坐标.txtfilenotopen!endl;exit(1);}while(f1m){x[i][j++]=m;if(j=3){i++;j=0;}}f1.close();cout左片各点像空间坐标坐标为:endl;for(i=0;i6;i++)for(j=0;j3;j++){if(j%3==0)coutendl;coutsetw(8)x[i][j];}coutendl;i=0;j=0;ifstreamf2(右片各点像空间坐标.txt);if(!f2){cerr右片像各点空间坐标.txtfilenotopen!endl;exit(1);}while(f2m){y[i][j++]=m;if(j=3){i++;j=0;}}f2.close();cout右片各点像空间坐标坐标为:endl;for(i=0;i6;i++)for(j=0;j3;j++){if(j%3==0)coutendl;coutsetw(8)y[i][j];}coutendl;b=0.000527;}voidmain(){doublea[5]={0},d[5]={0},R[3][3],N[2][6],A[6][5],AT[5][6],l[6][1],ATA[5][5],ATl[5][1],DG[5][1];intt=0;Input();do{t++;for(i=0;i5;i++)a[i]=a[i]+d[i];//计算旋转矩阵R[0][0]=cos(a[0])*cos(a[2])-sin(a[0])*sin(a[1])*sin(a[2]);R[0][1]=-cos(a[0])*sin(a[2])-sin(a[0])*sin(a[1])*cos(a[2]);R[0][2]=-sin(a[0])*cos(a[1]);R[1][0]=cos(a[1])*sin(a[2]);R[1][1]=cos(a[1])*cos(a[2]);R[1][2]=-sin(a[1]);R[2][0]=sin(a[0])*cos(a[2])+cos(a[0])*sin(a[1])*sin(a[2]);R[2][1]=-sin(a[0])*sin(a[2])+cos(a[0])*sin(a[1])*cos(a[2]);R[2][2]=cos(a[0])*cos(a[1]);//计算右片各点空间辅助坐标for(i=0;i6;i++)for(j=0;j3;j++)z[i][j]=R[j][0]*y[i][0]+R[j][1]*y[i][1]+R[j][2]*y[i][2];for(i=0;i6;i++){N[0][i]=(b*z[i][2]-b*a[4]*z[i][0])/(x[i][0]*z[i][2]-z[i][0]*x[i][2]);N[1][i]=(b*x[i][2]-b*a[4]*x[i][0])/(x[i][0]*z[i][2]-z[i][0]*x[i][2]);}for(i=0;i6;i++)l[i][0]=N[0][i]*x[i][1]-N[1][i]*z[i][1]-b*a[3];for(i=0;i6;i++){A[i][0]=-z[i][0]*z[i][1]/z[i][2]*N[1][i];A[i][1]=-(z[i][2]+z[i][1]*z[i][1]/z[i][2])*N[1][i];A[i][2]=z[i][0]*N[1][i];A[i][3]=b;A[i][4]=-z[i][1]*b/z[i][2];}Transpose(A,AT,6,5);Array_mul(AT,A,ATA,5,6,5);Inverse(ATA,5);Array_mul(AT,l,ATl,5,6,1);Array_mul(ATA,ATl,DG,5,5,1);for(i=0;i5;i++)for(j=0;j1;j++)d[i]=DG[i][j];}while((fabs(d[0])0.00003)||(fabs(d[1])0.00003)||(fabs(d[2])0.00003)||(fabs(d[3])0.00003)||(fabs(d[4])0.00003));cout迭代次数为:tendl;cout未知数的解为:endl;coutψ=a[0],ω=a[1],κ=a[2],;c