作业4空间后方交会Word文档下载推荐.docx
《作业4空间后方交会Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《作业4空间后方交会Word文档下载推荐.docx（21页珍藏版）》请在冰豆网上搜索。
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<
stdio.h>
stdlib.h>
math.h>
constdoublePRECISION=1e-5;
typedefdoubleDOUBLE[5];
intInputData（int&
Num,DOUBLE*&
Data,double&
m,double&
f）;
intResection（constint&
Num,constDOUBLE*&
Data,constdouble&
m,constdouble&
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!
printf（"
解算完毕...\n"
）;
do{
printf（"
计算结果保存于\"
结果.txt\"
文件中\n"
"
请选择操作（输入P打开结果数据，R打开原始数据，其它退出程序）：
"
fflush（stdin）;
//刷新输入流
charorder=getchar（）;
if（'
P'
==order||'
p'
==order）
{
system（"
结果.txt"
elseif（'
R'
r'
data.txt"
else
break;
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文件格式错误
**********************************************/
f）
doublex0,y0;
FILE*fp_input;
if（!
（fp_input=fopen（"
"
r"
）））
fscanf（fp_input,"
%d%lf"
&
Num,&
m）;
if（Num<
4）
return2;
%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<
Num;
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）
//单位换算成m
Data[i][0]/=1000.0;
Data[i][1]/=1000.0;
//如果m未知则归算其值
if（0==m）
for（i=0;
Num-1;
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;
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;
row;
for（intj=0;
j<
line;
j++）
doubletemp=0;
for（intk=0;
k<
com;
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*&
m0,
constdouble&
Xs,constdouble&
Ys,constdouble&
Zs,
Phi,constdouble&
Omega,
Kappa,constdouble*R）
FILE*fp_out;
（fp_out=fopen（"
w"
fprintf（fp_out,"
**************************************"
*********************************\n"
\n空间后方交会程序（C\\C++）\n测绘一班\n"
学号：
20080729\n姓名：
陈闻亚\n\n"
已知数据：
\n\n已知点数：
intnum;
doubletemp,x,y;
num,&
temp）;
%d\n"
num）;
摄影比例尺（0表示其值位置）：
%10.0lf\n"
temp）;
内方位元素（fx0y0）：
temp,&
x,&
y）;
%10lf\t%10lf\t%10lf\n"
temp,x,y）;
num;
doubletemp[5];
fscanf（fp_input,"
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]）;
计算结果如下：
\n\n外方位元素：
\n"
\tXs=%10lf\n"
Xs）;
\tYs=%10lf\n"
Ys）;
\tZs=%10lf\n"
Zs）;
\tPhi=%10lf\n"
Phi）;
\tOmega=%10lf\n"
Omega）;
\tKappa=%10lf\n\n"
Kappa）;
旋转矩阵：
3;
\t"
%10lf\t"
*（R+i*3+j））;
\n单位权中误差：
%10lf\n\n"
m0）;
协因数阵：
6;
%20lf\t"
*（Q+i*6+j））;
\n外方位元素精度：
m[i]）;
fclose（fp_out）;
Resection
计算
Num（点数），Data（数据），m,,f（焦距）,x0,y0
int，0成功，其它失败
m,
doubleXs=0,Ys=0,Zs=0;
inti,j;
//设置初始值
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）;
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）*
A[2*i+1][5]=-（xxx）;
//求矩阵A的转置矩阵AT
2*Num;
for（j=0;
*（AT+j*2*Num+i）=A[i][j];
//求ATA
MatrixMul（AT,6,&
A[0][0],6,2*Num,ATA）;
if（InverseMatrix（ATA,6））
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];
Kappa+=X
升级会员 