绝对定向作业报告

绝对定向作业报告1 作业任务------------------------------------------------------------------------------------ 32 作业原理--------------------------------------------------------------------------------------- 33 已知条件及数据-------------------------------------------------------------------- 34 作业过程--------------------------------------------------------------------------- 35 源程序----------------------------------------------------------------------------- 46 计算结果--------------------------------------------------------------------------- 10 7心得体会与建议----------------------------------------------------------------------------- 101 作业任务在上次作业（相对定向）的基础上，继续完成以下编程任务：1）计算出6个定向点的模型坐标、摄影测量坐标（航摄比例尺为1:37000）；2）采用如下绝对定向元素：l ＝1.00156, ψ＝0.0527, ω＝0.1426, κ=0.2478, X0=6385.067, Y0=1954.325, Z0=724.215，计算出各点的地面摄影测量坐标；3）根据地面摄影坐标和模型点的摄影坐标，编程实现绝对定向，计算出绝对定向元素和各点的地面摄影测量坐标。

2 作业原理解析法绝对定向：利用已知的地面控制点，从绝对定向的关系式出发，解求绝对定向元素。

实际上，绝对定向德主要工作室把模型点的摄影测量坐标变换为地面摄影测量坐标。

空间相似变换：个立体像对有12个外方位元素，经相对定向求得五个定向元素后，要恢复像对的绝对方位，还要求解7个绝对定向元素，包括模型的旋转、平移和缩放。

这种坐标变换前后图形的几何形状相似，称为空间相似变换。

设任一模型点的摄影测量坐标为（U,V,W ），对应的地面摄影测量坐标为（X,Y,Z ），它们之间的空间相似变换可以用绝对定向的基本关系式（下式）表示，即X Y Z ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦ = λ 123123123a a a b b b c c c ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦ U V W ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦ + s s s X Y Z ⎡⎤⎢⎥⎢⎥⎢⎥⎣⎦式中：λ为缩放系数；ai ，bi ，ci 为由角元素ψ, ω, κ的函数组成的方向余弦；Xs ，Ys ，Zs 为坐标原点的平移量。

解析绝对定向就是根据控制点的地面摄影测量坐标和对应的模型坐标（摄测坐标），解算出ψ, ω, κ, Xs , Ys , Zs 和λ共7个绝对定向参数，再用算得的7个参数，把待定点的摄影测量坐标换算为地面摄影测量坐标。

本次作业为先利用给定绝对定向元素求出坐标条件，再以求得的坐标条件为已知条件，又求绝对定向元素，与所给值作对比，并求出地面摄影测量坐标的改正值。

3已知条件及数据采用如下绝对定向元素：l ＝1.00156, ψ＝0.0527, ω＝0.1426, κ=0.2478, X0=6385.067, Y0=1954.325, Z0=724.215，计算出各点的地面摄影测量坐标后，又以求得的地面摄影测量坐标和根据上次作业解算的模型点摄影测量坐标为已知条件，代入绝对定向的基本关系式。

4作业过程4．1模型点坐标计算根据上次相对定向作业，正确求解出相对定向元素后，利用空间前方交会共识计算出模型点的坐标。

任一模型点坐标：111122Y 11m 1m b 2m X N X Y N Y N Y Z N Z =⎧⎪⎪=++⎨⎪=⎪⎩（）4．2 求模型点的摄影测量坐标：为了后续计算，把上步求得坐标平移到摄影测量坐标系中，同时放大模型比例尺，使之接近实地大小。

111122Y 111b 2U mN X V N Y N Y W mf mN Z =⎧⎪⎪=++⎨⎪=+⎪⎩m （）4．3根据所给绝对定向元素重新计算旋转矩阵R ，以及上述坐标共同代入绝对定向基本关系式，求解各点的地面摄影测量坐标（X,Y,Z ）。

4．4将所求得得各点的地面摄影测量坐标和模型点的摄影测量坐标作为已知条件，绝对定向元素视为未知，进行绝对定向。

4．5确定相对定向元素的初始值；ψ0 = ω0 = κ0 =0， Xs 0=Ys 0 =Zs 0 =0，λ0 =14．6根据确定的初始值（或新的近似值），计算出误差方程式的常数项。

4．7逐点组成误差方程并法化，逐点法化。

4．8解求法方程，得七个绝对定向元素的改正数。

4．8计算绝对定向元素的新值。

4．9判断绝对定向元素的改正数是否小于限值0.00003 rad ，如满足条件，则结束相对定向计算。

否则重复4.5~4.9。

4．10根据求得的绝对定向元素，将所有模型点的摄测坐标转换为地面摄测坐标。

5 源程序#include<iostream.h>#include<fstream.h>#include<stdlib.h>#include<iomanip.h>#include<math.h>const int N=7;int i,j,k;double b,x[6][3], y[6][3],z[6][3]={0};//x[6][3]和y[6][3]分别为左右片像点像空间坐标，z[6][3]为右片各点像空间辅助坐标//求转置矩阵template<typename T1,typename T2>void Transpose(T1*mat1,T2*mat2,int a,int b){for(i=0;i<a;i++)for(j=0;j<b;j++)mat2[j][i]=mat1[i][j];return;}//求矩阵的乘积template<typename T1,typename T2>void Array_mul(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c){ int i,j,k;for(i=0;i<a;i++){for(j=0;j<c;j++){result[i][j]=0;for(k=0;k<b;k++)result[i][j]+=mat1[i][k]*mat2[k][j];}}return;}//求逆矩阵void swap(double *a,double *b){double c; c= *a; *a= *b; *b= c;};Inverse(double A[N][N],int n){int i,j,k;double d;int JS[N],IS[N];for (k=0;k<n;k++){d=0;for (i=k;i<n;i++)for (j=k;j<n;j++){if (fabs(A[i][j])>d){d=fabs(A[i][j]);IS[k]=i;JS[k]=j;};};if (d+1.0==1.0) return 0;if (IS[k]!=k) for (j=0;j<n;j++) swap( &A[k][j], &A[IS[k]][j]);if (JS[k]!=k) for (i=0;i<n;i++) swap( &A[i][k], &A[i][JS[k]]);A[k][k]=1/A[k][k];for (j=0;j<n;j++) if (j!=k) A[k][j]=A[k][j]*A[k][k];for (i=0;i<n;i++) if (i!=k) for (j=0;j<n;j++) if (j!=k) A[i][j]=A[i][j]-A[i][k]*A[k][j];for (i=0;i<n;i++) if (i!=k) A[i][k]=-A[i][k]*A[k][k];};for (k=n-1;k>=0;k--){for (j=0;j<n;j++) if (JS[k]!=k) swap( &A[k][j], &A[JS[k]][j]);for (i=0;i<n;i++) if (IS[k]!=k) swap( &A[i][k], &A[i][IS[k]]);};return 1;}//原始数据导入void Input(){ double m;j=0;i=0;ifstream f1("左片各点像空间坐标.txt");if(!f1){cerr<<"左片各点像空间坐标.txt file not open!"<<endl;exit(1);}while(f1>>m){x[i][j++]=m;if(j>=3){i++;j=0;}}f1.close();cout<<"左片各点像空间坐标坐标为:"<<endl;for(i=0;i<6;i++)for(j=0;j<3;j++){if(j%3==0)cout<<endl;cout<<setw(8)<<x[i][j];}cout<<endl;i=0;j=0;ifstream f2("右片各点像空间坐标.txt");if(!f2){cerr<<"右片像各点空间坐标.txt file not open!"<<endl;exit(1);}while(f2>>m){y[i][j++]=m;if(j>=3){i++;j=0;}}f2.close();cout<<"右片各点像空间坐标坐标为:"<<endl;for(i=0;i<6;i++)for(j=0;j<3;j++){if(j%3==0)cout<<endl;cout<<setw(8)<<y[i][j];}cout<<endl;b=0.000527;}void main(){double a[5]={0},d[5]={0}, R[3][3], N[2][6],A[6][7]={0}, AT[7][6]={0}, l[6][1]={0},ATA[7][7]={0},ATl[7][1]={0},DG[7][1]={0};int t=0;Input();do{t++;for(i=0;i<5;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;i<6;i++)for(j=0;j<3;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;i<6;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;i<6;i++)l[i][0]=N[0][i]*x[i][1]-N[1][i]*z[i][1]-b*a[3];for(i=0;i<6;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;i<5;i++)d[i]=DG[i][0];}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<<"相对定向元素计算迭代次数为:"<<t<<endl;cout<<"未知数的解为:"<<endl;cout<<"ψ="<<a[0]<<","<<"ω="<<a[1]<<","<<"κ="<<a[2]<<",";cout<<"μ="<<a[3]<<","<<"γ="<<a[4]<<"."<<endl;ofstream f3("计算结果.txt");if(!f3){cerr<<"计算结果.txt file not open!"<<endl;exit(1);}f3<<"迭代次数为:"<<t<<endl;f3<<"相对定向元素的值为:"<<endl;f3<<"ψ="<<a[0]<<","<<"ω="<<a[1]<<","<<"κ="<<a[2]<<",";f3<<"bv="<<b*a[3]<<","<<"bw="<<b*a[4]<<"."<<endl;f3<<endl;f3<<endl;//计算出个像点的模型坐标mzb[6][3]、摄影测量坐标mszb[6][3]（m=37000），地面摄影测量坐标ds[6][3] //采用如下绝对定向元素：缩放系数k＝.00156, ψ＝.0527, ω＝.1426, κ=0.2478, X0=6385.067,Y0=1954.325, Z0=724.215，计算出各点的地面摄影测量坐标；double mzb[6][3]={0},mszb[6][3]={0},ds[6][3]={0},m=37000.0,f=0.024,k=1.00156,X0=6385.067, Y0=1954.325, Z0=724.215;a[0]=0.0527;a[1]=0.1426;a[2]=0.2478;for(i=0;i<6;i++){mzb[i][0]=N[0][i]*x[i][0];mzb[i][1]=0.5*(N[0][i]*x[i][1]+N[1][i]*y[i][1]+b*a[3]);mzb[i][2]=N[0][i]*x[i][2];}for(i=0;i<6;i++){mszb[i][0]=mzb[i][0]*m;mszb[i][1]=mzb[i][1]*m;mszb[i][2]=m*f+mzb[i][2]*m;}cout<<"6个像点的模型坐标分别为："<<endl;for(i=0;i<6;i++){cout<<mzb[i][0]<<" , "<<mzb[i][1]<<" , "<<mzb[i][2]<<"."<<endl;}cout<<"6个像点的摄影测量坐标分别为："<<endl;for(i=0;i<6;i++){cout<<mszb[i][0]<<" , "<<mszb[i][1]<<" , "<<mszb[i][2]<<"."<<endl;}//计算各点的地面摄影测量坐标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;i<6;i++){ds[i][0]=k*(R[0][0]*mszb[i][0]+R[0][1]*mszb[i][1]+R[0][2]*mszb[i][2])+X0;ds[i][1]=k*(R[1][0]*mszb[i][0]+R[1][1]*mszb[i][1]+R[1][2]*mszb[i][2])+Y0;ds[i][2]=k*(R[2][0]*mszb[i][0]+R[2][1]*mszb[i][1]+R[2][2]*mszb[i][2])+Z0;}cout<<"各点的地面摄影测量坐标为:"<<endl;for(i=0;i<6;i++){cout<<ds[i][0]<<" , "<<ds[i][1]<<" , "<<ds[i][2]<<"."<<endl;}f3<<"计算得各点的地面摄影测量坐标为:"<<endl;f3<<" X Y Z "<<endl;for(i=0;i<6;i++){f3<<ds[i][0]<<" , "<<ds[i][1]<<" , "<<ds[i][2]<<"."<<endl;}f3.close();//根据地面摄影坐标和模型点的摄影坐标，实现绝对定向，计算出绝对定向元素和各点的地面摄影测量坐标doublejX[7][1],dX[7],jA[18][7],jAT[7][18],jATA[7][7],jL[18][1],jATL[7][1],jG[7][1],jV[18][1],jAX[1 8][1];//重心化坐标double Xp=0,Yp=0,Zp=0,Xmp=0,Ymp=0,Zmp=0;for(i=0;i<6;i++){Xp=Xp+mszb[i][0];Yp=Yp+mszb[i][1];Zp=Zp+mszb[i][2];Xmp=Xmp+ds[i][0];Ymp=Ymp+ds[i][1];Zmp=Zmp+ds[i][2];}Xp=Xp/6;Yp=Yp/6;Zp=Zp/6;Xmp=Xmp/6;Ymp=Ymp/6;Zmp=Zmp/6;for(i=0;i<6;i++){mszb[i][0]=mszb[i][0]-Xp;mszb[i][1]=mszb[i][1]-Yp;mszb[i][2]=mszb[i][2]-Zp;ds[i][0]=ds[i][0]-Xmp;ds[i][1]=ds[i][1]-Ymp;ds[i][2]=ds[i][2]-Zmp;}t=0;do{t++;for(i=0;i<7;i++)jX[i][0]=jX[i][0]+dX[i];//计算旋转矩阵R[0][0]=cos(jX[4][0])*cos(jX[6][0])-sin(jX[4][0])*sin(jX[5][0])*sin(jX[6][0]);R[0][1]=-cos(jX[4][0])*sin(jX[6][0])-sin(jX[4][0])*sin(jX[5][0])*cos(jX[6][0]);R[0][2]=-sin(jX[4][0])*cos(jX[5][0]);R[1][0]=cos(jX[5][0])*sin(jX[6][0]);R[1][1]=cos(jX[5][0])*cos(jX[6][0]);R[1][2]=-sin(jX[5][0]);R[2][0]=sin(jX[4][0])*cos(jX[6][0])+cos(jX[4][0])*sin(jX[5][0])*sin(jX[6][0]);R[2][1]=-sin(jX[4][0])*sin(jX[6][0])+cos(jX[4][0])*sin(jX[5][0])*cos(jX[6][0]);R[2][2]=cos(jX[4][0])*cos(jX[5][0]);//计算常数项Lfor(i=0;i<6;i++){jL[i*3][0]=ds[i][0]-jX[3][0]*(R[0][0]*mszb[i][0]+R[0][1]*mszb[i][1]+R[0][2]*mszb[i][2])-jX[0 ][0];jL[i*3+1][0]=ds[i][1]-jX[3][0]*(R[1][0]*mszb[i][0]+R[1][1]*mszb[i][1]+R[1][2]*mszb[i][2])-jX [1][0];jL[i*3+2][0]=ds[i][2]-jX[3][0]*(R[2][0]*mszb[i][0]+R[2][1]*mszb[i][1]+R[2][2]*mszb[i][2])-jX [2][0];}//构造矩阵Afor(i=0;i<6;i++){jA[i*3][0]=1;jA[i*3][1]=0;jA[i*3][2]=0;jA[i*3][3]=mszb[i][0];jA[i*3][4]=-mszb[i][2];jA[i*3][ 5]=0;jA[i*3][6]=-mszb[i][1];jA[i*3+1][0]=0;jA[i*3+1][1]=1;jA[i*3+1][2]=0;jA[i*3+1][3]=mszb[i][1];jA[i*3+1][4]=0;jA[i*3+1][5]=-mszb[i][2];jA[i*3+1][6]=mszb[i][1];jA[i*3+2][0]=0;jA[i*3+2][1]=0;jA[i*3+2][2]=1;jA[i*3+2][3]=mszb[i][2];jA[i*3+2][4]=mszb[i][0] ;jA[i*3+2][5]=mszb[i][1];jA[i*3+2][6]=0;}Transpose(jA,jAT,18,7);Array_mul(jAT,jA,jATA,7,18,7);Inverse(jATA,7);Array_mul(jAT,jL,jATL,7,18,1);Array_mul(jATA,jATL,jG,7,7,1);for(i=0;i<7;i++)dX[i]=jG[i][0];}while((fabs(dX[4])>0.00003)||(fabs(dX[5])>0.00003)||(fabs(dX[6])>0.00003));//计算各点的地面摄影测量坐标改正数Array_mul(jA,jX,jAX,18,7,1);for(i=0;i<18;i++)jV[i][0]=jL[i][0]-jAX[i][0];//结果cout<<"迭代次数为:"<<t<<endl;cout<<"7个绝对定向元素的解为:"<<endl;cout<<"Xs="<<jX[0]<<","<<"Ys="<<jX[1]<<","<<"Zs="<<jX[2]<<"."<<endl;cout<<"λ="<<jX[3]<<"."<<endl;cout<<"ψ="<<jX[4]<<","<<"ω="<<jX[5]<<","<<"κ="<<jX[6]<<".";cout<<"各点的地面摄影测量坐标改正数为（从左到右分别为X,Y,Z）:"<<endl;for(i=0;i<6;i++){cout<<jV[i*3][0]<<" , "<<jV[i*3+1][0]<<" , "<<jV[i*3+2][0]<<endl;}ofstream f4("绝对定向元素计算结果.txt");if(!f4){cerr<<"绝对定向元素计算结果.txt file not open!"<<endl;exit(1);}f4<<"迭代次数为:"<<t<<endl;f4<<"绝对定向元素的值为:"<<endl;f4<<"Xs="<<jX[0]<<","<<"Ys="<<jX[1]<<","<<"Zs="<<jX[2]<<"."<<endl;f4<<"λ="<<jX[3]<<"."<<endl;f4<<"ψ="<<jX[4]<<","<<"ω="<<jX[5]<<","<<"κ="<<jX[6]<<"."<<endl;f4<<"各点的地面摄影测量坐标:"<<endl;for(i=0;i<6;i++){f4<<ds[i][0]<<"±"<<jV[i*3][0]<<" , "<<ds[i][1]<<"±"<<jV[i*3+1][0]<<" ,"<<ds[i][2]<<"±"<<jV[i*3+2][0]<<endl;}f4.close();}7 计算结果8 心得体会与建议这是本学期最后一次编程作业，此次作业总体来说是比较轻松的。