绝对定向作业报告

绝对定向作业报告

1 作业任务 ------------------------------------------------------------------------------------ 2 作业原理 --------------------------------------------------------------------------------------- 3 已知条件及数据-------------------------------------------------------------------- 4 作业过程 --------------------------------------------------------------------------- 5 源程序 -----------------------------------------------------------------------------

6 计算结果 --------------------------------------------------------------------------- 7心得体会与建议-----------------------------------------------------------------------------

3 3 3 3 4

10 10

1 作业任务

在上次作业（相对定向）的基础上，继续完成以下编程任务：

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a1a2a3Ub1b2b3 V + c1c2c3WXsYs Zs式中：λ为缩放系数；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模型点坐标计算

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

XmN1X11 Ym（N1Y1N2Y2bY）2ZmN1Z1

4．2 求模型点的摄影测量坐标：

为了后续计算，把上步求得坐标平移到摄影测量坐标系中，同时放大模型比例尺，使之接近实地大小。

UmN1X11 Vm（N1Y1N2Y2bY）2WmfmN1Z1

4．3根据所给绝对定向元素重新计算旋转矩阵R，以及上述坐标共同代入绝对定向基本关系式，求解各点的地面摄影测量坐标（X,Y,Z）。

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

4．5确定相对定向元素的初始值；

ψ0 = ω0 = κ0 =0， Xs0=Ys0 =Zs0 =0，λ0 =1

4．6根据确定的初始值（或新的近似值），计算出误差方程式的常数项。

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

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

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

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

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

5 源程序

#include #include #include #include #include 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]为右片各点像空间辅助坐标

//求转置矩阵

templatevoid Transpose(T1*mat1,T2*mat2,int a,int b) { }

//求矩阵的乘积

templatevoid Array_mul(T1*mat1,T2 * mat2,T2 * result,int a,int b,int c)

{ int i,j,k; }

//求逆矩阵

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;kfor(i=0;i{result[i][j]=0; }

for(k=0;kresult[i][j]+=mat1[i][k]*mat2[k][j];

for(i=0;ifor(j=0;jmat2[j][i]=mat1[i][j]; return;

for (i=k;id){ 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;jfor (j=0;jfor (i=0;ifor (k=n-1;k>=0;k--) {

for (j=0;jreturn 1;

} //原始数据导入 void Input()

{ double m;j=0;i=0;

ifstream f1(\"左片各点像空间坐标.txt\"); if(!f1) { }

while(f1>>m) { }

f1.close();

cout<<\"左片各点像空间坐标坐标为:\"<for(j=0;j<3;j++) x[i][j++]=m; if(j>=3) {i++;j=0;}

cerr<<\"左片各点像空间坐标.txt file not open!\"< }

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]);

{

if(j%3==0) cout<cout<cout<i=0;j=0;

ifstream f2(\"右片各点像空间坐标.txt\"); { }

while(f2>>m) { }

f2.close();

cout<<\"右片各点像空间坐标坐标为:\"<for(j=0;j<3;j++) {

if(j%3==0) cout<cout<cout<=3) {i++;j=0;}

cerr<<\"右片像各点空间坐标.txt file not open!\"<if(!f2)

for(i=0;i<6;i++)

for(i=0;i<6;i++)

Transpose(A,AT,6,5);

Inverse(ATA,5);

Array_mul(AT,A,ATA,5,6,5); Array_mul(AT,l,ATl,5,6,1); Array_mul(ATA,ATl,DG,5,5,1); for(i=0;i<5;i++)

}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<<\"相对定向元素计算迭代次数为:\"<cout<<\"ψ=\"<d[i]=DG[i][0]; {

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; }

l[i][0]=N[0][i]*x[i][1]-N[1][i]*z[i][1]-b*a[3];

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]); 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(j=0;j<3;j++)

for(i=0;i<6;i++)

z[i][j]=R[j][0]*y[i][0]+R[j][1]*y[i][1]+R[j][2]*y[i][2];

A[i][4]=-z[i][1]*b/z[i][2];

cout<<\"μ=\"<{cerr<<\"计算结果.txt file not open!\"<f3<<\"迭代次数为:\"<f3<<\"ψ=\"<//计算出个像点的模型坐标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个像点的模型坐标分别为：\"<{

cout<cout<<\"6个像点的摄影测量坐标分别为：\"<{

cout<//计算各点的地面摄影测量坐标

R[0][0]=cos(a[0])*cos(a[2])-sin(a[0])*sin(a[1])*sin(a[2]);

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; 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]);

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<<\"各点的地面摄影测量坐标为:\"<f3<<\"计算得各点的地面摄影测量坐标为:\"<f3<f3.close();

//根据地面摄影坐标和模型点的摄影坐标，实现绝对定向，计算出绝对定向元素和各点的地面摄影测量坐标

double

jX[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[18][1]; //重心化坐标

double Xp=0,Yp=0,Zp=0,Xmp=0,Ymp=0,Zmp=0; for(i=0;i<6;i++) {

Xmp=Xmp+ds[i][0];Ymp=Ymp+ds[i][1];Zmp=Zmp+ds[i][2]; } {

cout<Xp=Xp+mszb[i][0];Yp=Yp+mszb[i][1];Zp=Zp+mszb[i][2];

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]);

Xp=Xp/6;Yp=Yp/6;Zp=Zp/6;Xmp=Xmp/6;Ymp=Ymp/6;Zmp=Zmp/6; {

mszb[i][0]=mszb[i][0]-Xp;mszb[i][1]=mszb[i][1]-Yp;mszb[i][2]=mszb[i][2]-Zp; } t=0;

ds[i][0]=ds[i][0]-Xmp;ds[i][1]=ds[i][1]-Ymp;ds[i][2]=ds[i][2]-Zmp;

for(i=0;i<6;i++)

//计算常数项L for(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];

}

for(i=0;i<6;i++) {

//构造矩阵A

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++)

}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<<\"迭代次数为:\"<cout<<\"ψ=\"<cout<ofstream f4(\"绝对定向元素计算结果.txt\"); if(!f4)

{cerr<<\"绝对定向元素计算结果.txt file not open!\"<f4<<\"迭代次数为:\"<f4<<\"Xs=\"<f4<<\"ψ=\"<f4<dX[i]=jG[i][0];

}

}

f4.close(); }

7 计算结果

8 心得体会与建议

这是本学期最后一次编程作业，此次作业总体来说是比较轻松的。首先是在前几次作业的基础上，没有需要太多的变化，其次是我们自己对这类问题有了一定的掌握，解决起来也比较熟练。通过本次作业，使我真正地弄明白了相对定向、绝对定向的意义以及它们之间的区别和联系，心中不会再有忐忑的模糊感。总体来说，本次作业基本顺利，其中在编程方面遇到的一点小问题就是：矩阵求逆的子函数的两次调用。不过在与同学的沟通和讨论下，解决了这个问题。