摄影测量实习报告单片空间后方交会.docx

摄影测量实习报告单片空间后方交会.docx

《摄影测量实习报告单片空间后方交会.docx》由会员分享，可在线阅读，更多相关《摄影测量实习报告单片空间后方交会.docx（10页珍藏版）》请在冰豆网上搜索。

摄影测量实习报告单片空间后方交会

摄影测量实习报告

实习内容:

单片空间后方交会编程

实习者:

李友兵

学号:

**********

指导老师:

张金平老师

实习时间:

2011.05.30——2011.06.03

一、实习目的与任务

此次摄影测量实习主要是要自主编程实现单像空间后方交会，通过已知的内方位元素和控制点像点坐标和地面坐标求解六个外方位元素，在此过程中深入理解单像空间后方交会的原理和对编程的熟悉和理解（我用的是C语言编程），和对时间的合理运用，对知识的综合运用，培养理论的实际运用能力，任务是在一个星期内自主完成。

二、单片空间后方交会理论基础

单像空间后方交会：

是通过以像点平面坐标为观测值，以控制点坐标为已知值，利用共线条件方程和最小二乘原理，运用间接平差方法，通过迭代求解6个外方位元素。

程序设计的思路与流程

三、程序设计的思路与流程（编程框架和步骤）

1、根据内方位元素和已知的数据将控制点的框标坐标转换为像平面坐标系坐标：

x=x′-x0，y=y′-y0

2、确定未知数的初始值：

角元素初始值κ0=ω0=Φ0=0；线元素初始值Zs0=H=mf，Xs0=（X1+X2+X3+X4）/4,Ys0=（Y1+Y2+Y3+Y4）/4

3、利用角元素初始值计算方向余弦值组成旋转矩阵R

a1=cosΦ*cosκ-sinΦ*sinω*sinκ,a2=-cosΦ*sinκ-sinΦ*sinω*cosκ，a3=-sinΦ*cosω

b1=cosω*sinκ,b2=cosω*cosκ,b3=-sinω

c1=sinΦ*cosκ+cosΦ*sinω*sinκ,c2=-sinΦ*sinκ+cosΦ*sinω*cosκ,c3=cosΦ*cosω

4、逐点计算控制点的像点坐标的近似值，共线方程为:

x=（-f*（a1*（X-Xs）+b1*（Y-Ys）+c1*（Z-Zs）））/（a3*（X-Xs）+b3*（Y-Ys）+c3*（Z-Zs））

y=（-f*（a2*（X-Xs）+b2*（Y-Ys）+c2*（Z-Zs）））/（a3*（X-Xs）+b3*（Y-Ys）+c3*（Z-Zs））

（x）=（-f*（a1*（X-Xs）+b1*（Y-Ys）+c1*（Z-Zs）））/（a3*（X-Xs）+b3*（Y-Ys）+c3*（Z-Zs））

（y）=（-f*（a2*（X-Xs）+b2*（Y-Ys）+c2*（Z-Zs）））/（a3*（X-Xs）+b3*（Y-Ys）+c3*（Z-Zs））

5、组成误差方程：

Vx=a11*dXs+a12*dYs+a13*dZs+a14*dΦ+a15*dω+a16*dκ+（x）-x

Vy=a21*dXs+a22*dYs+a23*dZs+a24*dΦ+a25*dω+a26*dκ+（y）-y

系数求解（是共线方程分别外方位元素求导，是共线方程线性化的系数）：

变量代换

A=a1*（X-Xs）+b1*（Y-Ys）+c1*（Z-Zs）

B=a2*（X-Xs）+b2*（Y-Ys）+c2*（Z-Zs）

C=a3*（X-Xs）+b3*（Y-Ys）+c3*（Z-Zs）

a11=（a1*f+a3*x）/C,a12=（b1*f+b3*x）/C,a13=（c1*f+c3*f）/C,a14=y*sinω-（（x*（x*cosκ-y*sinκ））/f+f*cosκ）cosω,a15=-f*sinκ-（x*（x*sinκ+y*cosκ）/f）,a16=y

a21=（a2*f+a3*x）/C,a22=（b2*f+b3*x）/C,a23=（c2*f+c3*x）/C,a24=-xsinω-（（x*（x*cosκ-y*sinκ））/f-f*sinκ）cosω,a25=-f*cosκ-（y*（x*sinκ+y*cosκ）/f）,a26=-x

误差方程的常系数（是像点坐标观测值与计算的近似值的差值）：

lx=x-（x）,ly=y-（y）

6、组成法方程，解求外方位元素改正数X=（ATA）-1ATL（A为误差方程的系数矩阵，L为误差方程的常系数矩阵，通过步骤5求得，此处先求ATA再求矩阵的逆矩阵，解得的改正数加上相应的近似值得到外方位元素新的近似值）

7、检查计算是否收敛：

将求得的外方位元素改正数与规定的限差比较，大于限差继续迭代，小于限差则终止。

四、各子函数详细设计的关键技术参数

子函数（输入函数、Input，矩阵求积Matrixmultiply，计算函数Resection，矩阵转置Matrixtranspose，矩阵求逆Matrixinverse，输出函数Output）主要用了Matrixtranspose，矩阵的行变列，列变行,参数为要转置的矩阵，转置后的矩阵，要转置矩阵的行列数，Matrixinverse,求矩阵的代数余子式，参数有要求逆的矩阵和，逆矩阵的行数，Matrixmultiply，一矩阵的行乘以二矩阵的列，参数为一矩阵，二矩阵，所求矩阵，一的行，一的列，二的列

五、像片外方位元素解算结

六、实习体会

从周一开始进行整个实习框架进行构建，编写了实习的思路既步骤；然后是梳理框架，要用到的子函数，子函数参数，怎么编写，有不懂的东西网上参考资料；再是自己动手编程，发现编程要求很细腻，出不得一点差错，程序编写出来后，有很多的错误，要求逐一进行改正，修改；最后才得以运行。

通过本次实习，我深刻理解了单片空间后方交会原理，进一步熟悉理解了C语言编程，认识到了时间的搭配的重要性和参考资料的必要性，当遇到难题是要迎难而上，达到突破，当完成是能感到一丝丝欣慰和成就感。

附录:

源代码

#include"stdio.h"

#include"math.h"

#include"Matrixmultiply.c"

#include"Matrixtranspose.c"

#include"Matrixinverse.c"

voidmain（）

{

inti,j,k,f=0;

doublex0=0.00018,y0=0.00026,fk=0.15324;//内方位元素

doublem=40000;//估算比例尺

doubleR[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];

doubleXs=0.0,Ys=0.0,Zs=0.0,Q=0.0,W=0.0,K=0.0;

doubleX,Y,Z,L[8][1],A[8][6];

doubleB[4][5]={-0.08615,-0.06899,36589.41,25273.32,2195.17,-0.05340,0.08221,37631.08,31324.51,728.69,-0.01478,-0.07663,39100.97,24934.98,2386.80,0.01046,0.06443,40426.54,30319.81,757.31};

for（i=0;i<4;i++）

{

Xs=Xs+B[i][2];

Ys=Ys+B[i][3];

Zs=Zs+B[i][4];

}

Xs=Xs/4;Ys=Ys/4;Zs=m*fk;//求得外方位线元素的初始值

do//迭代计算

{

f++;//迭代次数

//组成旋转矩阵

R[0][0]=cos（Q）*cos（K）-sin（Q）*sin（W）*sin（K）;

R[0][1]=-cos（Q）*sin（K）-sin（Q）*sin（W）*cos（K）;

R[0][2]=-sin（Q）*cos（W）;

R[1][0]=cos（W）*sin（K）;

R[1][1]=cos（W）*cos（K）;

R[1][2]=-sin（W）;

R[2][0]=sin（Q）*cos（K）+cos（Q）*sin（W）*sin（K）;

R[2][1]=-sin（Q）*sin（K）+cos（Q）*sin（W）*cos（K）;

R[2][2]=cos（Q）*cos（W）;

//计算系数阵和常数项

for（i=0,k=0,j=0;i<=3;i++,k++,j++）

{

X=R[0][0]*（B[i][2]-Xs）+R[1][0]*（B[i][3]-Ys）+R[2][0]*（B[i][4]-Zs）;

Y=R[0][1]*（B[i][2]-Xs）+R[1][1]*（B[i][3]-Ys）+R[2][1]*（B[i][4]-Zs）;

Z=R[0][2]*（B[i][2]-Xs）+R[1][2]*（B[i][3]-Ys）+R[2][2]*（B[i][4]-Zs）;//为了计算简单而进行变量代换

L[j][0]=B[i][0]-（x0-fk*X/Z）;

L[j+1][0]=B[i][1]-（y0-fk*Y/Z）;

j++;

A[k][0]=（R[0][0]*fk+R[0][2]*（B[i][0]-x0））/Z;

A[k][1]=（R[1][0]*fk+R[1][2]*（B[i][0]-x0））/Z;

A[k][2]=（R[2][0]*fk+R[2][2]*（B[i][0]-x0））/Z;

A[k][3]=（B[i][1]-y0）*sin（W）-（（B[i][0]-x0）*（（B[i][0]-x0）*cos（K）-（B[i][1]-y0）*sin（K））/fk+fk*cos（K））*cos（W）;

A[k][4]=-fk*sin（K）-（B[i][0]-x0）*（（B[i][0]-x0）*sin（K）+（B[i][1]-y0）*cos（K））/fk;

A[k][5]=B[i][1]-y0;

A[k+1][0]=（R[0][1]*fk+R[0][2]*（B[i][1]-y0））/Z;

A[k+1][1]=（R[1][1]*fk+R[1][2]*（B[i][1]-y0））/Z;

A[k+1][2]=（R[2][1]*fk+R[2][2]*（B[i][1]-y0））/Z;

A[k+1][3]=-（B[i][0]-x0）*sin（W）-（（B[i][1]-y0）*（（B[i][0]-x0）*cos（K）-（B[i][1]-y0）*sin（K））/fk-fk*sin（K））*cos（W）;

A[k+1][4]=-fk*cos（K）-（B[i][1]-y0）*（（B[i][0]-x0）*sin（K）+（B[i][1]-y0）*cos（K））/fk;

A[k+1][5]=-（B[i][0]-x0）;

k++;

}

Matrixtranspose（A,AT,8,6）;//此为转置函数的调用,求AT

Matrixmultiply（AT,A,ATA,6,8,6）;//此为矩阵相乘函数的调用，求ATA

Matrixinverse（ATA,6）;//此为求矩阵逆函数的调用，求ATA的逆

Matrixmultiply（AT,L,ATL,6,8,1）;//此为矩阵相乘函数的调用，求ATL

Matrixmultiply（ATA,ATL,XG,6,6,1）;//此为矩阵相乘函数的调用，求外方位元素改正数

Xs=Xs+XG[0][0];Ys=Ys+XG[1][0];Zs=Zs+XG[2][0];

Q=Q+XG[3][0];W=W+XG[4][0];K=K+XG[5][0];//初始值加外方位元素改正数进行迭代

}while（XG[3][0]>=0.00000001||XG[4][0]>=0.00000001||XG[5][0]>=0.00000001）;

//当限差满足要求时要再一次进行旋转矩阵的求解

R[0][0]=cos（Q）*cos（K）-sin（Q）*sin（W）*sin（K）;

R[0][1]=-cos（Q）*sin（K）-sin（Q）*sin（W）*cos（K）;

R[0][2]=-sin（Q）*cos（W）;

R[1][0]=cos（W）*sin（K）;

R[1][1]=cos（W）*cos（K）;

R[1][2]=-sin（W）;

R[2][0]=sin（Q）*cos（K）+cos（Q）*sin（W）*sin（K）;

R[2][1]=-sin（Q）*sin（K）+cos（Q）*sin（W）*cos（K）;

R[2][2]=cos（Q）*cos（W）;

printf（"迭代次数:

%d",f）;

//屏幕输出误差方程系数阵、常数项、改正数

printf（"\n\n误差方程系数矩阵A为:

\n\n"）;

for（i=0;i<6;i++）

{

for（j=0;j<6;j++）

printf（"%13.5e",A[i][j]）;

printf（"\n"）;

}

printf（"\n常数项L为:

\n\n"）;

for（i=0;i<8;i++）

{

for（j=0;j<1;j++）

printf（"%13.5e",L[i][j]）;

printf（"\n"）;

}

printf（"\n改正数XG为:

\n\n"）;

for（i=0;i<6;i++）

{

for（j=0;j<1;j++）

printf（"%13.5e",XG[i][j]）;

printf（"\n"）;

}

printf（"\n相片的外方位元素为:

\n\n"）;

printf（"Xs=%13.7e,Ys=%13.7e,Zs=%13.7e\n\n",Xs,Ys,Zs）;

printf（"Q=%13.7e,W=%13.7e,K=%13.7e\n",Q,W,K）;

printf（"\n旋转矩阵R为:

\n\n"）;

for（i=0;i<3;i++）

{

for（j=0;j<3;j++）

printf（"%13.5e",R[i][j]）;

printf（"\n"）;

}

}

//子函数

#include

#include

#include

intMatrixinverse（a,n）

intn;

doublea[];

{int*is,*js,i,j,k,l,u,v;

doubled,p;

is=malloc（n*sizeof（int））;

js=malloc（n*sizeof（int））;

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）

{free（is）;free（js）;printf（"err**notinv\n"）;

return（0）;

}

if（is[k]!

=k）

for（j=0;j<=n-1;j++）

{u=k*n+j;v=is[k]*n+j;

p=a[u];a[u]=a[v];a[v]=p;

}

if（js[k]!

=k）

for（i=0;i<=n-1;i++）

{u=i*n+k;v=i*n+js[k];

p=a[u];a[u]=a[v];a[v]=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;v=js[k]*n+j;

p=a[u];a[u]=a[v];a[v]=p;

}

if（is[k]!

=k）

for（i=0;i<=n-1;i++）

{u=i*n+k;v=i*n+is[k];

p=a[u];a[u]=a[v];a[v]=p;

}

}

free（is）;free（js）;

return

（1）;

}

//子函数

voidMatrixmultiply（a,b,c,m,n,k）

intm,n,k;

doublea[],b[],c[];

{

inti,j,l,u;

for（i=0;i

for（j=0;j

{

u=i*k+j;c[u]=0.0;

for（l=0;l

c[u]+=a[i*n+l]*b[l*k+j];

}

return;

}

//子函数

voidMatrixtranspose（a,b,m,n）

intm,n;

doublea[],b[];

{

inti,j,u;

for（i=0;i

{

for（j=0;j

{

u=j*m+i;b[u]=0.0;

b[j*m+i]=a[i*n+j];

}

}

return;

}