单张影像空间后方交会程序.docx
- 文档编号:6549401
- 上传时间:2023-05-10
- 格式:DOCX
- 页数:17
- 大小:299.03KB
单张影像空间后方交会程序.docx
《单张影像空间后方交会程序.docx》由会员分享,可在线阅读,更多相关《单张影像空间后方交会程序.docx(17页珍藏版)》请在冰点文库上搜索。
单张影像空间后方交会程序
摄影测量实习报告
-----单张影像空间后方交会程序
实习时间2013.5.20-2013.5.24
学生班级10测绘
(1)班
学生姓名房新明
学生学号**********
所在院系矿业工程学院
指导老师张会战邵亚琴
一.实习目的
1.深入理解单张影像空间后方交会的原理,体会在有多余观测情况下,用最小二乘平差方法编程实现解求影像外方位元素的过程。
2.利用VisualC++编写一个完整的单张影像空间后方交会程序,通过对提供的试验数据进行计算,输出像片的外方位元素并进行评定精度。
3.通过编写程序实现单张影像空间后方交会计算,掌握非线性方程线性化的过程、相应数据读入与存储的方法以及迭代计算的特点,巩固各类基础课程及计算机课程的学习内容,培养上机调试程序的动手能力,通过对实验结果的分析,增强综合运用所学知识解决专业实际问题的能力。
4、掌握空间后方交会的定义和实现算法
(1)定义:
空间后方交会是以单幅影像为基础,从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发,根据共线条件方程,解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,φ,ω,κ。
(2)算法:
由于每一对像方和物方共轭点可列出2个方程,因此若有3个已知地面坐标的控制点,则可列出6个方程,解求6个外方位元素的改正数△Xs,△Ys,△Zs,△φ,△ω,△κ。
实际应用中为了提高解算精度,常有多余观测方程,通常是在影像的四个角上选取4个或均匀地选择更多的地面控制点,因而要用最小二乘平差方法进行计算。
5、了解空间后方交会的基本过程
(1)获取已知数据。
从摄影资料中查取影像比例尺1/m,平均摄影距离(航空摄影的航高)、内方位元素x0,y0,f;获取控制点的空间坐标Xt,Yt,Zt。
(2)量测控制点的像点坐标并进行必要的影像坐标系统误差改正,得到像点坐标。
(3)确定未知数的初始值。
单像空间后方交会必须给出待定参数的初始值,在竖直航空摄影且地面控制点大体对称分布的情况下,Xs0和Ys0为均值,Zs0为航高,φ、ω、κ的初值都设为0。
或者κ的初值可在航迹图上找出或根据控制点坐标通过坐标正反变换求出。
(4)计算旋转矩阵R。
利用角元素近似值计算方向余弦值,组成R阵。
(5)逐点计算像点坐标的近似值。
利用未知数的近似值按共线条件式计算控制点像点坐标的近似值(x),(y)。
(6)逐点计算误差方程式的系数和常数项,组成误差方程式。
(7)计算法方程的系数矩阵ATA与常数项ATL,组成法方程式。
(8)解求外方位元素。
根据法方程,解求外方位元素改正数,并与相应的近似值求和,得到外方位元素新的近似值。
(9)检查计算是否收敛。
将所求得的外方位元素的改正数与规定的限差比较,通常对φ,ω,κ的改正数△φ,△ω,△κ给予限差,通常为0.1′,当3个改正数均小于0.1′时,迭代结束。
否则用新的近似值重复(4)~(8)步骤的计算,直到满足要求为止。
二、实习要求
1、认真复习单张影像空间后方交会的原理及有关内容;
2、编写程序,计算左片号23和右片号24的外方位元素并进行精度评定;
3、验证数据部分外方位元素和旋转矩阵R已经给出,请验证自己编写的程序是否正确;
三、实习环境
1、硬件环境:
windows操作系统
2、软件环境:
VC++
四、实习原理
1、以单幅影像为基础,从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发,根据共线条件方程,求解该影象在航空摄影时刻的像片外方位元素Xs,Ys,zs,∮,ω,κ共线条件方程如下:
x-x0=-f*[a1(X-Xs)+b1(Y-Ys)+c1(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)]
y-y0=-f*[a2(X-Xs)+b2(Y-Ys)+c2(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)]
其中:
x,y为像点的像平面坐标;Xs,Ys,f为影像的外方位元素;
Xs,Ys,Zs为摄站点的物方空间坐标;
X,Y,Z为物方点的物方空间坐标;
旋转矩阵为R;
由于此共线条件方程是非线性方程,先对其进行线性化,像点观测值一般视为等权,即P=I;
矩阵形式:
V=AX-L,P=I;
通过间接平差,为提高精度,增加多余观测方程,根据最小二乘平差原理,可计算出外方位元素的改正数。
经过迭代计算,每次迭代用未知数的近似值与上次迭代计算的改正数之和作为新的近似值,重复计算,求出新的改正数,这样反复趋近,直到改正数小于某个限值为止。
2、精度评定
Mi=m0*√Qii其中m0=±√【VV】/(2N-6)
3、公式推导
六、程序过程框图:
七.实习数据
1、模拟像片一对:
左片23号右片24号
2、相片比例尺:
1/30000
3、航摄机主距:
f=150mm
4、每张相片有四个控制点
5、各片像点坐标及地面坐标
片号
点号
像点坐标(mm)
地面坐标(m)
x
y
X
Y
Z
23
1
-91.596
-74.859
100000.0
137500.0
11.00
3
-94.230
81.446
100000.0
142500.00
36.00
7
95.207
-75.512
106000.0
137500.00
42.00
9
96.797
83.077
106000.0
142500.00
56.00
24
4
-102.695
-79.618
103000.0
137500.00
90.00
6
-99.904
81.754
103000.0
142500.00
31.00
10
86.890
-77.540
109000.0
137500.00
7.00
12
88.904
76.257
109000.0
142500.00
5.00
八.验证数据
已知四对点的影像坐标和地面坐标及航摄仪内方位元素f=153.42mm,X0=Y0=0
影像坐标
地面坐标
x(mm)
Y(mm)
X(mm)
Y(mm)
Z(mm)
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
3、结算结果
Xs=39795.450.997710.067530.00399
Ys=27476.46R=-0.067530.99772-0.00221
Zs=7572.69-0.004120.001840.99990
九.实现程序
输入文件形式如下:
C++源程序如下:
#include
#include
#include
#include
#include
usingnamespacestd;
constintn=6;
voidinverse(doublec[n][n]);
template
template
template
template
intmain()
{
ofstreamoutFile;
cout.precision(5);
doublex0=0.0,y0=0.0;doublefk=0.15324;//内方位元素
doublem=39689;//估算比例尺
doubleB[4][5]={0.0},R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];
input(B,4,5);//从文件中读取控制点的影像坐标和地面坐标,存入数组B
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];
//确定未知数的出始值
for(inti=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=Zs/4+m*fk;
intf=0;
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(inti=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++;
}
transpose(A,AT,6,8);
multi(AT,A,ATA,6,8,6);
inverse(ATA);
multi(AT,L,ATL,6,8,1);
multi(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]>=6.0/206265.0||XG[4][0]>=6.0/206265.0||XG[5][0]>=6.0/206265.0);
cout<<"迭代次数为:
"< //精度评定 doubleAXG[8][1],V[8][1],VT[1][8],VTV[1][1],m0,D[6][6]; multi(A,XG,AXG,8,6,1); for(i=0;i<8;i++)//计算改正数 V[i][0]=AXG[i][0]-L[i][0]; transpose(V,VT,1,8); multi(VT,V,VTV,1,8,1); m0=VTV[0][0]/2; for(i=0;i<6;i++) for(intj=0;j<6;j++) D[i][j]=m0*ATA[i][j]; //屏幕输出误差方程系数阵、常数项、改正数 output(A,"误差方程系数阵A为: ",8,6); output(L,"常数项L为: ",8,1); output(XG,"改正数为: ",6,1); outFile.open("aim.txt",ios: : app);//打开并添加aim.txt文件 outFile.precision(10); //以文件的形式输出像片外方位元素、旋转矩阵、方差阵 outFile<<"一、像片的外方位元素为: "< outFile< < outFile< "< "旁向倾角为: "< "< outFile< "< for(i=0;i<3;i++) {for(intj=0;j<3;j++) outFile< outFile< } outFile< outFile< "< outFile.precision(5); for(i=0;i<6;i++) {for(intj=0;j<6;j++) outFile< outFile< } outFile.close(); return0; } template {inti,j; for(i=0;i for(j=0;j mat2[j][i]=mat1[i][j]; return; } template {inti,j,k; for(i=0;i {for(j=0;j {result[i][j]=0; for(k=0;k result[i][j]+=mat1[i][k]*mat2[k][j]; } } return; } template {ifstreaminFile; inFile.open("控制点坐标.txt"); while(! inFile.eof()) {for(inti=0;i for(intj=0;j inFile>>mat[i][j]; } inFile.close(); return; } template
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 单张 影像 空间 后方 交会 程序