作业4空间后方交会.docx
《作业4空间后方交会.docx》由会员分享，可在线阅读，更多相关《作业4空间后方交会.docx（21页珍藏版）》请在冰豆网上搜索。
作业4空间后方交会
作业报告
空间后方交会
专业：
测绘工程
班级：
2008级
（1）班
姓名：
陈闻亚
指导教师：
陈强
2010年4月16日
1作业任务------------------------------------------------------------------------------------3
2作业思想---------------------------------------------------------------------------------------3
3作业条件及数据--------------------------------------------------------------------3
4作业过程---------------------------------------------------------------------------3
5源程序-----------------------------------------------------------------------------4
6计算结果---------------------------------------------------------------------------17
7心得体会与建议-----------------------------------------------------------------------------17
1作业任务
计算近似垂直摄影情况下后方交会解。
即利用摄影测量空间后方交会的方法，获取相片的6个外方位元素。
限差为0.1。
2作业思想
利用摄影测量空间后方交会的方法求解。
该方法的基本思想是利用至少三个一直地面控制点的坐标A（XA，YA，ZA）、B（XB，YB，ZB）C（XC，YC，ZC），与其影像上对应的三个像点的影像坐标a（xa，ya）、b（xb，yb）、c（xc，yc），根据共线方程，反求该相片的外方位元素XS、YS、ZS、φ、ω、κ。
3作业条件及数据
已知摄影机主距f=153.24mm，四对点的像点坐标与相应的地面坐标列入下表：
表1
点号
像点坐标
地面坐标
x（mm）
y（mm）
X（m）
Y（m）
Z（m）
1
-86.15
-68.99
36589.41
25273.32
2195.17
2
-53.40
-82.21
37631.08
31324.51.
728.69
3
14.78
76.63
39100.97
24934.98
2386.50
4
10.46
64.43
40426.54
30319.81
757.31
4作业过程
4．1获取已知数据
相片比例尺1/m=1：
10000，内方位元素f=153.24mm，x0，y0；获取控制点的地面测量坐标Xt、Yt、Zt。
4．2量测控制点的像点坐标：
本次作业中为已知。
见表1。
4．3确定未知数的初始值：
在近似垂直摄影情况下，胶原素的初始值为0，即φ0=ω0=κ0=0；线元素中，ZS0=H=mf=1532.4m，XS0、YS0的取值可用四个控制点坐标的平均值，即：
XS0=
=38437.00
YS0=
=89106.62
4．4计算旋转矩阵R：
利用胶原素的近似值计算方向余弦值，组成R阵。
4．5逐点计算像点坐标的近似值：
利用未知数的近似值按共线方程式计算控制点像点坐标的近似值（x）（y）。
4．6组成误差方程：
逐点计算误差方程式的系数和常数项。
4．7组成法方程式：
计算法方程的系数矩阵ATA与常数项ATL。
4．8求解外方位元素：
根据法方程，由式X=（AtA）-1ATL解求外方位元素改正数，并与相应的近似值求和，得到外方位元素新的近似值。
4．9求解外方位元素：
将求得的外方位元素的改正数与规定的限差（0.1）比较，小于限差则计算终止，否则用新的近似值重复第4.4至4.8步骤的计算，知道满足要求为止。
5源程序
#include
#include
#include
constdoublePRECISION=1e-5;
typedefdoubleDOUBLE[5];
intInputData（int&Num,DOUBLE*&Data,double&m,double&f）;
intResection（constint&Num,constDOUBLE*&Data,constdouble&m,constdouble&f）;
intInverseMatrix（double*matrix,constint&row）;
intmain（intargc,char*argv[]）
{
DOUBLE*Data=NULL;
intNum;
doublef（0）,m（0）;
if（InputData（Num,Data,m,f））
{
if（Data!
=NULL）
{
delete[]Data;
}
return1;
}
if（Resection（Num,Data,m,f））
{
if（Data!
=NULL）
{
delete[]Data;
}
return1;
}
if（Data!
=NULL）
{
delete[]Data;
}
printf（"解算完毕...\n"）;
do{
printf（"计算结果保存于\"结果.txt\"文件中\n"
"请选择操作（输入P打开结果数据，R打开原始数据，其它退出程序）：
"）;
fflush（stdin）;//刷新输入流
charorder=getchar（）;
if（'P'==order||'p'==order）
{
system（"结果.txt"）;
}
elseif（'R'==order||'r'==order）
{
system（"data.txt"）;
}
else
break;
system（"cls"）;
}while
（1）;
system（"PAUSE"）;
return0;
}
/**********************************************
*函数名：
InputData
*函数介绍：
从文件（data.txt）中读取数据,
*文件格式如下：
*点数m（未知写作0）
*内方位元素（fx0y0）
*编号xyXYZ
*实例：
40
153.2400
1-86.15-68.9936589.4125273.322195.17
2-53.4082.2137631.0831324.51728.69
3-14.78-76.6339100.9724934.982386.50
410.4664.4340426.5430319.81757.31
*参数：
（in/out）Num（点数），
*（in/out）Data（存放数据）,m,f,x0,y0
*返回值：
int，0成功，1文件打开失败，2控制点个
*数不足，3文件格式错误
**********************************************/
intInputData（int&Num,DOUBLE*&Data,double&m,double&f）
{
doublex0,y0;
FILE*fp_input;
if（!
（fp_input=fopen（"data.txt","r"）））
{
return1;
}
fscanf（fp_input,"%d%lf",&Num,&m）;
if（Num<4）
{
return2;
}
fscanf（fp_input,"%lf%lf%lf",&f,&x0,&y0）;
f/=1000;
if（m<0||f<0）
{
return3;
}
Data=newDOUBLE[Num];
double*temp=newdouble[Num-1];
doublescale=0;
inti;
for（i=0;i{
//读取数据，忽略编号
if（fscanf（fp_input,"%*d%lf%lf%lf%lf%lf",
&Data[i][0],&Data[i][1],&Data[i][2],
&Data[i][3],&Data[i][4]）!
=5）
{
return3;
}
//单位换算成m
Data[i][0]/=1000.0;
Data[i][1]/=1000.0;
}
//如果m未知则归算其值
if（0==m）
{
for（i=0;i{
temp[i]=（Data[i][2]-Data[i+1][2]）/（Data[i][0]-Data[i+1][0]）+
（Data[i][3]-Data[i+1][3]）/（Data[i][1]-Data[i+1][1]）;
scale+=temp[i]/2.0;
}
m=scale/（Num-1）;
}
fclose（fp_input）;
delete[]temp;
return0;
}
/**********************************************
*函数名：
MatrixMul
*函数介绍：
求两个矩阵的积,
*参数：
Jz1（第一个矩阵）,row（第一个矩阵行数）,
*Jz2（第二个矩阵）,row（第二个矩阵列数）,com（第一个
*矩阵列数）,（out）JgJz（存放结果矩阵）
*返回值：
void
**********************************************/
voidMatrixMul（double*Jz1,constint&row,double*Jz2,
constint&line,constint&com,double*JgJz）
{
for（inti=0;i{
for（intj=0;j{
doubletemp=0;
for（intk=0;k{
temp+=*（Jz1+i*com+k）*（*（Jz2+k*line+j））;
}
*（JgJz+i*line+j）=temp;
}
}
}
/**********************************************
*函数名：
OutPut
*函数介绍：
向结果.txt文件输出数据
*参数：
Q协因数阵，m精度，m0单位权中误差，6个外
*方位元素，旋转矩阵
*返回值：
int，0成功，1失败
**********************************************/
intOutPut（constdouble*&Q,constdouble*&m,constdouble&m0,
constdouble&Xs,constdouble&Ys,constdouble&Zs,
constdouble&Phi,constdouble&Omega,
constdouble&Kappa,constdouble*R）
{
FILE*fp_out;
if（!
（fp_out=fopen（"结果.txt","w"）））
{
return1;
}
FILE*fp_input;
if（!
（fp_input=fopen（"data.txt","r"）））
{
return1;
}
fprintf（fp_out,"**************************************"
"**************************************"
"**************************************"
"*********************************\n"）;
fprintf（fp_out,"\n空间后方交会程序（C\\C++）\n测绘一班\n"
"学号：
20080729\n姓名：
陈闻亚\n\n"）;
fprintf（fp_out,"**************************************"
"**************************************"
"**************************************"
"*********************************\n"）;
fprintf（fp_out,"已知数据：
\n\n已知点数：
"）;
intnum;doubletemp,x,y;
fscanf（fp_input,"%d%lf",&num,&temp）;
fprintf（fp_out,"%d\n",num）;
fprintf（fp_out,"摄影比例尺（0表示其值位置）：
"）;
fprintf（fp_out,"%10.0lf\n",temp）;
fprintf（fp_out,"内方位元素（fx0y0）：
"）;
fscanf（fp_input,"%lf%lf%lf",&temp,&x,&y）;
fprintf（fp_out,"%10lf\t%10lf\t%10lf\n",temp,x,y）;
for（inti=0;i{
doubletemp[5];
fscanf（fp_input,"%*d%lf%lf%lf%lf%lf",
&temp[0],&temp[1],&temp[2],&temp[3],&temp[4]）;
fprintf（fp_out,"%3d\t%10lf\t%10lf\t%10lf\t%10lf\t%10lf\n",
i+1,temp[0],temp[1],temp[2],temp[3],temp[4]）;
}
fclose（fp_input）;
fprintf（fp_out,"**************************************"
"**************************************"
"**************************************"
"*********************************\n"）;
fprintf（fp_out,"计算结果如下：
\n\n外方位元素：
\n"）;
fprintf（fp_out,"\tXs=%10lf\n",Xs）;
fprintf（fp_out,"\tYs=%10lf\n",Ys）;
fprintf（fp_out,"\tZs=%10lf\n",Zs）;
fprintf（fp_out,"\tPhi=%10lf\n",Phi）;
fprintf（fp_out,"\tOmega=%10lf\n",Omega）;
fprintf（fp_out,"\tKappa=%10lf\n\n",Kappa）;
fprintf（fp_out,"旋转矩阵：
\n"）;
for（i=0;i<3;i++）
{
fprintf（fp_out,"\t"）;
for（intj=0;j<3;j++）
{
fprintf（fp_out,"%10lf\t",*（R+i*3+j））;
}
fprintf（fp_out,"\n"）;
}
fprintf（fp_out,"\n单位权中误差：
%10lf\n\n",m0）;
fprintf（fp_out,"协因数阵：
\n"）;
for（i=0;i<6;i++）
{
fprintf（fp_out,"\t"）;
for（intj=0;j<6;j++）
{
fprintf（fp_out,"%20lf\t",*（Q+i*6+j））;
}
fprintf（fp_out,"\n"）;
}
fprintf（fp_out,"\n外方位元素精度：
"）;
for（i=0;i<6;i++）
{
fprintf（fp_out,"%10lf\t",m[i]）;
}
fprintf（fp_out,"\n"）;
fprintf（fp_out,"**************************************"
"**************************************"
"**************************************"
"*********************************\n"）;
fclose（fp_out）;
return0;
}
/**********************************************
*函数名：
Resection
*函数介绍：
计算
*参数：
Num（点数），Data（数据），m,,f（焦距）,x0,y0
*返回值：
int，0成功，其它失败
**********************************************/
intResection（constint&Num,constDOUBLE*&Data,constdouble&m,
constdouble&f）
{
doubleXs=0,Ys=0,Zs=0;
inti,j;
//设置初始值
for（i=0;i{
Xs+=Data[i][2];
Ys+=Data[i][3];
}
Xs/=Num;
Ys/=Num;
Zs=m*f;
doublePhi（0）,Omega（0）,Kappa（0）;
doubleR[3][3]={0.0};
double*L=newdouble[2*Num];
typedefdoubleDouble6[6];
Double6*A=newDouble6[2*Num];
double*AT=newdouble[2*Num*6];
double*ATA=newdouble[6*6];
double*ATL=newdouble[6];
double*Xg=newdouble[6];
//迭代计算
do
{
//旋转矩阵
R[0][0]=cos（Phi）*cos（Kappa）-sin（Phi）*sin（Omega）*sin（Kappa）;
R[0][1]=-cos（Phi）*sin（Kappa）-sin（Phi）*sin（Omega）*cos（Kappa）;
R[0][2]=-sin（Phi）*cos（Omega）;
R[1][0]=cos（Omega）*sin（Kappa）;
R[1][1]=cos（Omega）*cos（Kappa）;
R[1][2]=-sin（Omega）;
R[2][0]=sin（Phi）*cos（Kappa）+cos（Phi）*sin（Omega）*sin（Kappa）;
R[2][1]=-sin（Phi）*sin（Kappa）+cos（Phi）*sin（Omega）*cos（Kappa）;
R[2][2]=cos（Phi）*cos（Omega）;
for（i=0;i{
doubleX=R[0][0]*（Data[i][2]-Xs）+R[1][0]*（Data[i][3]-Ys）+
R[2][0]*（Data[i][4]-Zs）;
doubleY=R[0][1]*（Data[i][2]-Xs）+R[1][1]*（Data[i][3]-Ys）+
R[2][1]*（Data[i][4]-Zs）;
doubleZ=R[0][2]*（Data[i][2]-Xs）+R[1][2]*（Data[i][3]-Ys）+
R[2][2]*（Data[i][4]-Zs）;
doublexxx,yyy;
xxx=-f*X/Z;
yyy=-f*Y/Z;
//常数项
L[2*i]=Data[i][0]-（-f*X/Z）;
L[2*i+1]=Data[i][1]-（-f*Y/Z）;
A[2*i][0]=（R[0][0]*f+R[0][2]*（xxx））/Z;
A[2*i][1]=（R[1][0]*f+R[1][2]*（xxx））/Z;
A[2*i][2]=（R[2][0]*f+R[2][2]*（xxx））/Z;
A[2*i][3]=（yyy）*sin（Omega）-（（（xxx）/f）*
（（xxx）*cos（Kappa）-（yyy）*sin（Kappa））+
f*cos（Kappa））*cos（Omega）;
A[2*i][4]=-f*sin（Kappa）-（（xxx）/f）*（（xxx）*
sin（Kappa）+（yyy）*cos（Kappa））;
A[2*i][5]=（yyy）;
A[2*i+1][0]=（R[0][1]*f+R[0][2]*（yyy））/Z;
A[2*i+1][1]=（R[1][1]*f+R[1][2]*（yyy））/Z;
A[2*i+1][2]=（R[2][1]*f+R[2][2]*（yyy））/Z;
A[2*i+1][3]=-（xxx）*sin（Omega）-（（（yyy）/f）*
（（xxx）*cos（Kappa）-（yyy）*sin（Kappa））-
f*sin（Kappa））*cos（Omega）;
A[2*i+1][4]=-f*cos（Kappa）-（（yyy）/f）*（（xxx）*
sin（Kappa）+（yyy）*cos（Kappa））;
A[2*i+1][5]=-（xxx）;
}
//求矩阵A的转置矩阵AT
for（i=0;i<2*Num;i++）
{
for（j=0;j<6;j++）
{
*（AT+j*2*Num+i）=A[i][j];
}
}
//求ATA
MatrixMul（AT,6,&A[0][0],6,2*Num,ATA）;
if（InverseMatrix（ATA,6））
return1;
MatrixMul（AT,6,L,1,2*Num,ATL）;
MatrixMul（ATA,6,ATL,1,6,Xg）;
Xs+=Xg[0];
Ys+=Xg[1];
Zs+=Xg[2];
Phi+=Xg[3];
Omega+=Xg[4];
Kapp
升级会员 