httpreadpudncomdownloads74sourcecodewindowscsharp269472wdqkjhfjh空间后方交会实习报告docWord文档下载推荐.docx
- 文档编号:6797975
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:15
- 大小:32.70KB
httpreadpudncomdownloads74sourcecodewindowscsharp269472wdqkjhfjh空间后方交会实习报告docWord文档下载推荐.docx
《httpreadpudncomdownloads74sourcecodewindowscsharp269472wdqkjhfjh空间后方交会实习报告docWord文档下载推荐.docx》由会员分享,可在线阅读,更多相关《httpreadpudncomdownloads74sourcecodewindowscsharp269472wdqkjhfjh空间后方交会实习报告docWord文档下载推荐.docx(15页珍藏版)》请在冰点文库上搜索。
-76.63
39100.97
24934.98
2386.50
4
10.46
64.43
40426.54
30319.81
757.31
5实习过程:
5.1学习单张像片空间后方交会的基本理论,掌握其基本思想。
如果我们知道每幅影像的6个外方位元素,就能确定被摄物体与航摄影像的关系。
而单像空间后方交会就是用于测定像片的外方位元素的,它的基本思想是:
以单幅影像为基础,从影像所覆盖的地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发,根据共线方程,解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,p,w,k.由于空间后方交会所采用的数学模型共线方程是非线性函数,为了便于外方位元素的解求,首先将其线性化。
5.2在纸上绘出空间后方交会的计算机程序框图。
为了能够在宏观上指导我们编写程序,我们需要在草稿纸上绘出程序框图。
框图如下:
输入原始数据
↓
归算像点坐标x,y
↓
计算和确定初值Xs0,Ys0,Zs0,p0,w0,k0
︱组成旋转矩阵R
︱↓
︱计算(x),(y)和lx,ly
迭↓
迭逐点组成误差方程式并法化
次↓
数否所有点完否?
小↓完
于解法方程,求未知数改正数
限↓
差计算改正后的外方位元素
否否↓
?
未知数改正数<
限差否?
否︱↓是
︱整理并输出计算结果
输出中间结果和出错信息↓
︱正常结束
非正常结束
6.按照程序框图编写程序。
程序代码如下:
//#include"
stdafx.h"
#include<
math.h>
stdlib.h>
stdio.h>
#defineq206265
//矩阵相乘
voidmult(double*m1,double*m2,double*result,intm,intp,intn)
{
inti,j,k;
for(i=0;
i<
m;
i++)
for(j=0;
j<
n;
j++)
{
result[i*n+j]=0.0;
for(k=0;
k<
p;
k++)
result[i*n+j]+=m1[i*p+k]*m2[j+k*n];
}
return;
}
//矩阵求逆
intinvers_matrix(double*m1,intn)
int*is,*js;
inti,j,k,l,u,v;
doubletemp,max_v;
is=(int*)malloc(n*sizeof(int));
js=(int*)malloc(n*sizeof(int));
if(is==NULL||js==NULL)
{
printf("
outofmemory!
\n"
);
return(0);
}
for(k=0;
max_v=0.0;
for(i=k;
for(j=k;
temp=fabs(m1[i*n+j]);
if(temp>
max_v)
{
max_v=temp;
is[k]=i;
js[k]=j;
}
if(max_v==0.0)
free(is);
free(js);
printf("
inversisnotavailble!
return(0);
if(is[k]!
=k)
for(j=0;
u=k*n+j;
v=is[k]*n+j;
temp=m1[u];
m1[u]=m1[v];
m1[v]=temp;
if(js[k]!
for(i=0;
{
u=i*n+k;
v=i*n+js[k];
temp=m1[u];
}
l=k*n+k;
m1[l]=1.0/m1[l];
for(j=0;
if(j!
{
u=k*n+j;
m1[u]*=m1[l];
}
for(i=0;
if(i!
for(j=0;
if(j!
{
u=i*n+j;
m1[u]-=m1[i*n+k]*m1[k*n+j];
}
for(i=0;
if(i!
{
u=i*n+k;
m1[u]*=-m1[l];
}
}
for(k=n-1;
k>
=0;
k--)
{
if(js[k]!
u=k*n+j;
v=js[k]*n+j;
temp=m1[u];
for(i=0;
u=i*n+k;
v=i*n+is[k];
free(is);
free(js);
return
(1);
//结果输出格式
voidstarpic()
{printf("
\t\t\t\t1.程序过程\t\t\t\n"
\t\t\t\t\t\t\t\t\n"
\t\t\t\t2.结果\t\t\t\n"
\t\t\t\t3.退出\t\t\t\n"
voidmain()
{//赋初值
inti,j,k=0,flag,num=0,d=5,Num;
doubleM,x0,y0,f,x1,y1;
doubleXs=0,Ys=0,Zs=0,aq=0,wq=0,kq=0,v=0,m;
doubleR[3][3],A[8][6],At[6][8],X[6][1],L[8][1],Ai[6][6],Ait[6][8],precision[6],V[8][1];
doublea1,a2,a3,b1,b2,b3,c1,c2,c3,X0,Y0,Z0;
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.50,0.01046,0.06443,40426.54,30319.81,757.31};
starpic();
loop:
\n\t\t请输入你要进行操作的选项:
"
scanf("
%d"
&
Num);
while(Num!
=1&
&
Num!
=2&
=3)
printf("
\t\t你的输入有误!
!
请输入你要进行操作的选项:
if(Num==1)
{
输入比例尺分母M:
\n"
%lf"
M);
输入像主点坐标x0,y0:
scanf("
%lf%lf"
x0,&
y0);
输入摄影机主距f:
f);
4;
请按下面的顺序输入您的第%d组数据:
像点坐标x,y对应的地面控制点坐标X,Y,Z\n"
i+1);
for(j=0;
5;
B[i][j]);
if(Num==2)
M=50000,x0=0,y0=0,f=0.15324;
if(Num==3)
\n\t\t\t*!
exit(-1);
Xs+=B[i][2];
Ys+=B[i][3];
//控制点坐标取平均值
Xs=Xs/4;
Ys=Ys/4;
Zs=M*f;
do
{//开始循环,9个方向余弦
a1=R[0][0]=cos(aq)*cos(kq)-sin(aq)*sin(wq)*sin(kq);
a2=R[0][1]=-cos(aq)*sin(kq)-sin(aq)*sin(wq)*cos(kq);
a3=R[0][2]=-sin(aq)*cos(wq);
b1=R[1][0]=cos(wq)*sin(kq);
b2=R[1][1]=cos(wq)*cos(kq);
b3=R[1][2]=-sin(wq);
c1=R[2][0]=sin(aq)*cos(kq+cos(aq)*sin(wq)*sin(kq));
c2=R[2][1]=-sin(aq)*sin(kq)+cos(aq)*sin(wq)*cos(kq);
c3=R[2][2]=cos(aq)*cos(wq);
{//由共线方程求像点坐标近似值
X0=a1*(B[i][2]-Xs)+b1*(B[i][3]-Ys)+c1*(B[i][4]-Zs);
Y0=a2*(B[i][2]-Xs)+b2*(B[i][3]-Ys)+c2*(B[i][4]-Zs);
Z0=a3*(B[i][2]-Xs)+b3*(B[i][3]-Ys)+c3*(B[i][4]-Zs);
//A矩阵中的元素计算
A[2*i][0]=(a1*f+a3*(B[i][0]-x0))/Z0;
A[2*i][1]=(b1*f+b3*(B[i][0]-x0))/Z0;
A[2*i][2]=(c1*f+c3*(B[i][0]-x0))/Z0;
A[2*i][3]=(B[i][1]-y0)*sin(wq)-(B[i][0]-x0)*(B[i][0]-x0)*cos(kq)*cos(wq)/f-f*cos(kq)*cos(wq)+(B[i][0]-x0)*(B[i][1]-y0)*sin(kq)*cos(kq)/f;
A[2*i][4]=-f*sin(kq)-(B[i][0]-x0)*(B[i][0]-x0)*sin(kq)/f-(B[i][0]-x0)*(B[i][1]-y0)*cos(kq)/f;
A[2*i][5]=B[i][1]-y0;
A[2*i+1][0]=(a2*f+a3*(B[i][1]-y0))/Z0;
A[2*i+1][1]=(b2*f+b3*(B[i][1]-y0))/Z0;
A[2*i+1][2]=(c2*f+c3*(B[i][1]-y0))/Z0;
A[2*i+1][3]=-(B[i][0]-x0)*sin(wq)-(B[i][1]-y0)*(B[i][0]-x0)*cos(kq)*cos(wq)/f+f*sin(kq)*cos(wq)+(B[i][1]-y0)*(B[i][1]-y0)*sin(kq)*cos(wq)/f;
A[2*i+1][4]=-f*cos(kq)-(B[i][1]-y0)*(B[i][0]-x0)*sin(kq)/f-(B[i][1]-y0)*(B[i][1]-y0)*cos(kq)/f;
A[2*i+1][5]=x0-B[i][0];
x1=-f*X0/Z0;
y1=-f*Y0/Z0;
//得到L矩阵
L[2*i][0]=B[i][0]-x1;
L[2*i+1][0]=B[i][1]-y1;
//转置
8;
6;
At[j][i]=A[i][j];
mult(*At,*A,*Ai,6,8,6);
if(invers_matrix(*Ai,6)==NULL)
exit(-1);
mult(*Ai,*At,*Ait,6,6,8);
mult(*Ait,*L,*X,6,8,1);
mult(*A,*X,*V,8,6,1);
Xs+=X[0][0];
Ys+=X[1][0];
Zs+=X[2][0];
aq+=X[3][0];
wq+=X[4][0];
kq+=X[5][0];
num++;
//迭带次数增加
flag=0;
for(i=3;
{//迭代控制条件
if((X[i][0]>
(0.1*3.1415/180/60))&
(num<
=d))
flag=1;
while(flag==1);
V[i][0]-=L[i][0];
2;
B[i][j]+=V[k][0];
k++;
v+=(V[i][0]*V[i][0]);
m=sqrt(v/(2*4-6));
//未知数中误差
precision[i]=m*sqrt(Ai[i][i]);
//单位权中误差
printf("
\n\n\n单片空间后方交会结果\n\n"
\t迭代次数=%d\n\n"
num);
\t像主点坐标\tx0=%.0f,y0=%.0f\n\n"
x0,y0);
\t比例尺分母\tM=%.0f\n\n"
M);
\t摄影机主距(单位:
m)\tf=%f\n\n"
f);
\t像点坐标和对应的地面控制点坐标(单位:
m):
\t第%d组点"
\t"
{if(j==0||j==1)
%f"
B[i][j]);
else
printf("
%.2f"
\t单位权中误差为(单位:
微米):
\tm=%.2f\n\n"
m*1e6);
\t该航摄相片的外方位元素为(单位:
米、弧度):
\n\tXs\tYs\t\tZs\taq\twq\t\tkq\t\t\n\t%.2f%.2f%.2f%f%f%f\n\n"
Xs,Ys,Zs,aq,wq,kq);
\t外方位元素对应的精度为:
\t\n\tXs\tYs\tZs\taq\twq\tkq\t\n"
precision[i]);
\n\n\n"
\t\t\t输出R矩阵:
\n\n"
3;
{for(j=0;
\t%.5f"
R[i][j]);
\n操作成功!
继续吗?
(1.需要2.不需要):
b2);
if(b2==1)
{starpic();
gotoloop;
\n\n\t\n"
\n\t\t\t\t谢谢使用!
7程序结果显示(略)
8实习心得与总结:
这次实习我很早就在准备,首先是反复地温习单片空间后方交会的理论知识和原理,即利用一定数量的地面控制点,根据共线条件方程,反求像片的外方位元素。
在本次实习中,我们要求用四个点地面控制点,用最小二乘法平差计算。
编程的过程中碰到问题也是很多的,首先是请教编程熟练的同学关于矩阵求逆,转置的编法,然后是调试和大量错误的修改。
最后又因为精度不够而使用别的方法。
用逐点法化代替最先的大矩阵一起运算等等。
受益很多。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- httpreadpudncomdownloads74sourcecodewindowscsharp269472wdqkjhfjh 空间 后方 交会 实习 报告 doc
链接地址:https://www.bingdoc.com/p-6797975.html