摄影测量实习报告.docx
《摄影测量实习报告.docx》由会员分享，可在线阅读，更多相关《摄影测量实习报告.docx（17页珍藏版）》请在冰豆网上搜索。
摄影测量实习报告
《解析摄影测量课程编程实习》
实习报告
学院:
遥感信息工程学院
班级:
1504班
学号:
2015302590021
姓名:
汤舒畅
实习地点:
201机房
指导教师:
李欣
2017年10月27日
1实习目的
用程序设计语言（VisualC++或者C语言），配置OpenCV库后，编写一个完整的单像空间后方交会程序，通过对提供的试验数据进行计算，输出像片在航空摄影时刻的外方位元素，并评定精度。
本实验的目的在于让学生深入理解单像空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。
通过上机调试程序加强动手能力的培养，并通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。
2实习原理与实习流程
2.1实习原理
如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。
因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。
可采取的方法有：
利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会，简称单像空间后方交会。
单像空间后方交会的基本思想是：
以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，经过一系列的步骤解求该影像在航空摄影时刻的外方位元素Xs，Ys，Zs，t，w，k。
2.2实习流程
（1）获取已知数据。
从摄影资料中查取影像比例尺1/m，平均摄影距离（航空摄影的航高、内方位元素x0，y0，f；获取控制点的空间坐标Xt，Yt，Zt。
在本次实习中，已知数据等均由资料查询或老师提供。
（2）量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。
在本次实习中，像点坐标可以从资料中查询得到。
得到像点坐标后输入原始数据。
（3）确定未知数的初始值。
单像空间后方交会必须给出待定参数的初始值，在竖直航空摄影且地面控制点大体对称分布的情况下，可按如下方法确定初始值：
式中：
m为摄影比例尺分母，n为控制点个数；
（4）计算旋转矩阵R。
利用角元素的近似值按相关计算方向余弦a1,a2,a3,b1,b2,b3,c1,c2,c3，组成R阵。
（5）逐点计算像点坐标的近似值。
利用未知数的近似值按共线条件方程计算控制点像点坐标的近似值（x）、（y）。
（6）按下式逐点计算误差方程式的系数和常数项，组成误差方程。
（7）计算法方程的系数矩阵ATA与常数项ATL，组成法方程；
（8）解求外方位元素。
根据法方程，解求外方位元素的改正数，并与相应的近似值求和，得到外方位元素新的近似值。
（9）检查计算是否收敛。
将所求得的外方位元素的改正数与规定的限差比较，通常对t、ω、κ、Xs、Ys、Zs的改正数Δt，Δω，Δκ，ΔXs，ΔYs，ΔZs给予限差，当改正数小于限差时，迭代结束。
否则用新的近似值重复（4）——（8）步骤计算，直到满足要求为止。
（10）空间后方交会的精度估计：
按上述方法所求得的影像外方位元素的精度可以通过法方程式中未知数的系数矩阵的逆阵（ATA）-1来解求，此时视像点坐标为等精度不相关观测值。
因为ATA）-1中第i个主对角线上的元素Qii就是法方程式中第i个未知数的权倒数，若单位权中误差为m0，则第i个未知数的中误差为：
mi=
0
当参加空间后方交会的控制点有n个时，则单位权中误差可按下式计算：
3实习数据、实习结果及其分析
3.1实习数据
本次实习的实验数据如下：
1）已知航摄仪的内方位元素：
fk=153.24mm，x0=y0=0.0mm，摄影比例尺为1:
15000；
2）4个地面控制点的地面坐标及其对应像点的像片坐标如下：
表一控制点地面坐标与像点坐标
点号
像片坐标（mm）
地面点坐标（m）
x
y
X
Y
Z
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.2实习结果及其分析
本次实习的实验结果如下：
图一实习结果展示
根据程序运行后得到的实验结果，对比着教材中提供的初值，可以得到以下结论：
（1）Xs，Ys，Zs得到的实验结果与教材提供的结果差值分别为0.02m，0.08m，0.06m，均存在差值但是差距较小。
而对于旋转矩阵，得到的实验结果与教材提供的理论值近似一致，因此这里近似的认为实验结果正确。
（2）总循环次数较少，仅有四次。
说明初值的设立较为理想，使得循环次数较少，减少计算量。
（3）影像坐标和地面点坐标在循环之后几乎不发生改变，说明在程序中的误差矩阵V结果很小，近似的接近于0。
而从单位权中误差的结果也能说明这种情况。
（4）对代码接下来的改进部分：
确定外方位元素结果值与理论值的差值来源，消除差值。
4实习小结
本次实习虽然难度整体不大，但是在实习的过程中我有了许多的收获，现总结如下：
（1）本次实习让我对单像空间后方交会这一摄影测量学科中非常重要的知识点有了更深的理解，通过影像中若干（≥3）个控制点的地面坐标和对应的像点坐标的量测值出发，我们可以利用单像空间后方交会来求解出影像拍摄时刻的六个外方位元素。
而具体的实现步骤，仅仅凭借着老师上课的讲解，不能完全的理解。
更需要课下实习，通过实现代码来更深刻的理解。
（2）OpenCV库可以比较好的保存，处理，输出C++程序中的矩阵类数据，而OpenCV的自带函数也可以很好地对矩阵数据进行运算。
通过OpenCV可以避免减少在矩阵处理函数编写方面许多不必要的精力，从而将精力更多的放在程序编写和精度的控制上。
（3）不同的初始值可能会导致不同的最终结果。
我原本对于Zs0，我在没有参考教材时，采用的是四个地面控制点Z坐标的均值，作为初始值。
此时得到的外方位元素Xs，Ys，Zs的结果，与教材中参考结果值分别相差1cm，8cm，15cm，而修改初始值，采用教材上讲解的方式后，最好得到的外方位元素Xs，Ys，Zs的结果与教材中参考值分别相差2cm，6cm，6cm。
精度整体上得到了提升。
但是为什么会出现这一种情况，我暂时还没有思路。
（4）本次实习也暴露出了我时间分配不合理的缺点，本次实习在上个月就已经布置了，但是我在接近一个月的时间中没有及时的分配时间，导致最后一个星期一直在赶作业，即使到交纸质版的前一页，我仍然在修改实习报告。
这一次也是对我敲响了一个警钟，我下次需要注意这一个问题。
5附件
本次实习的主要代码与讲解如下：
5.1头文件，公共变量和main（）函数部分
Main函数是本次实习的主函数，主要负责
（1）整个实习过程中体系框架的建立
（2）初始值的确定（3）最后程序结果的输出。
为了简化矩阵运算，本次实习配置了OpenCV，利用OpenCV中的Mat类数据结构和相关的矩阵操作函数，大大简化了程序。
本次实习中设置了一些公共变量，这些变量会在大部分的函数中作为已知值和改变值出现，因此在函数申明前，设立一些公共变量，以简化程序内容。
Main函数的最后部分将公共变量中的指针结构delete，释放出被占用的空间。
具体代码如下：
#include"stdafx.h"
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
usingnamespacestd;
usingnamespacecv;
float*x;
float*y;
float*X;
float*Y;
float*Z;
floatfk=153.24;
floatm=15000.0;
floatx00=0;
floaty00=0;
floato=0;//单位权中误差
float*mo;//各参数的精度
inttotal_number;
intmain（）
{
floata=0;//注意航向倾角为角度制
floatw=0;//注意旁侧倾角为角度制，需要转化为弧度制
floatk=0;//注意相片倾角为角度制，需要转化为弧度制
floatXs=（36589.41+37631.08+39100.97+40426.54）/4;
floatYs=（25273.32+31324.51+24934.98+30319.81）/4;
//floatZs=15000*fk+（2195.17+728.69+2389.50+757.50）/4;
floatZs=m*fk;
//单位权中误差精度赋值
mo=newfloat[6];
for（inti=0;i<6;i++）
{
mo[i]=0;
}
//储存坐标值
if（SaveTheCorrdinate（））
{
//计算矩阵差
intoutput=0;
output=CalculateV（&a,&w,&k,&Xs,&Ys,&Zs）;
intnumber=0;
while（output!
=1&&number<200）
{
number++;
output=CalculateV（&a,&w,&k,&Xs,&Ys,&Zs）;
}
FILE*fp=fopen（"answer.txt","w"）;
fprintf（fp,"影像坐标与地面点坐标如下：
\n"）;
fprintf（fp,"i\tx（mm）\ty（mm）\tX（m）\t\tY（m）\t\tZ（m）\n"）;
for（inti=1;i<=total_number;i++）
{
fprintf（fp,"%d\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n",i,x[i-1],y[i-1],X[i-1],Y[i-1],Z[i-1]）;
}
fprintf（fp,"单位权中误差如下：
\n"）;
fprintf（fp,"o=%.5f\n",o）;
fprintf（fp,"外方位元素及其精度如下：
\n"）;
fprintf（fp,"Xs=%.2f\tYs=%.2f\tZs=%.2f\n",Xs,Ys,Zs）;
fprintf（fp,"mXs=%.5f\tmYs=%.5f\tmZs=%.5f\n",mo[0],mo[1],mo[2]）;
fprintf（fp,"a=%.5f\tw=%.5f\tk=%.5f\n",a,w,k）;
fprintf（fp,"ma=%.5f\tmw=%.5f\tmk=%.5f\n",mo[3],mo[4],mo[5]）;
fprintf（fp,"总循环次数为：
%d",number+1）;
fclose（fp）;
waitKey（0）;
}
delete[]x;
delete[]y;
delete[]X;
delete[]Y;
delete[]Z;
delete[]mo;
return0;
}
5.2SaveTheCorrdinate（）部分
该部分的主要任务为从txt文本中读取原始数据，并将公共变量初始化，将原始数据存储在公共变量中。
boolSaveTheCorrdinate（）
{
//读取原始数据行数
ifstreamcorrdinatefile;
corrdinatefile.open（"coordinate.txt"）;
if（!
corrdinatefile）
{
printf（"Error!
Can'topenthecorrdinate!
\n"）;
returnfalse;
}
stringsline;
total_number=0;
while（getline（corrdinatefile,sline））
{
total_number++;
}
corrdinatefile.close（）;
x=newfloat[total_number];
y=newfloat[total_number];
X=newfloat[total_number];
Y=newfloat[total_number];
Z=newfloat[total_number];
//一维数组存储数据
corrdinatefile.open（"coordinate.txt"）;
total_number=0;
while（getline（corrdinatefile,sline））
{
istringstreamsin（sline）;
sin>>x[total_number]>>y[total_number]>>X[total_number]>>Y[total_number]>>Z[total_number];
total_number++;
}
corrdinatefile.close（）;
returntrue;
}
5.3CalculateV（）部分
这一部分是代码中主体计算部分，主要的任务是计算单像空间后方交会。
具体的部分如下：
（1）根据上一步循环结果，计算当前旋转矩阵
（2）根据旋转矩阵结果，计算当前贡献条件方程结果
（3）形成误差方程式并法化
（4）求解外方位元素改正数
（5）检查迭代是否收敛
具体代码如下：
intCalculateV（float*a,float*w,float*k,float*Xs,float*Ys,float*Zs）
{
//根据当前初值计算旋转矩阵
cv:
:
MatR=Mat（3,3,CV_32F）;
R.at（0,0）=1*cos（*a）*cos（*k）-sin（*a）*sin（*w）*sin（*k）;//a1
R.at（0,1）=-1*cos（*a）*sin（*k）-sin（*a）*sin（*w）*cos（*k）;//a2
R.at（0,2）=-1*sin（*a）*cos（*w）;//a3
R.at（1,0）=1*cos（*w）*sin（*k）;//b1
R.at（1,1）=1*cos（*w）*cos（*k）;//b2
R.at（1,2）=-1*sin（*w）;//b3
R.at（2,0）=1*sin（*a）*cos（*k）+cos（*a）*sin（*w）*sin（*k）;//c1
R.at（2,1）=-1*sin（*a）*sin（*k）+cos（*a）*sin（*w）*cos（*k）;//c2
R.at（2,2）=1*cos（*a）*cos（*w）;//c3
for（inti=0;i<3;i++）
{
for（intj=0;j<3;j++）
{
R_total[i][j]=R.at（i,j）;
}
}
//计算A矩阵之前先计算X1Y1Z1的值（代替X一八，Y一八，Z一八）
float*X1=newfloat[total_number];
for（inti=0;i{
X1[i]=0;
X1[i]=X1[i]+R.at（0,0）*（X[i]-（*Xs））;
X1[i]=X1[i]+R.at（1,0）*（Y[i]-（*Ys））;
X1[i]=X1[i]+R.at（2,0）*（Z[i]-（*Zs））;
}
float*Y1=newfloat[total_number];
for（inti=0;i{
Y1[i]=0;
Y1[i]=Y1[i]+R.at（0,1）*（X[i]-（*Xs））;
Y1[i]=Y1[i]+R.at（1,1）*（Y[i]-（*Ys））;
Y1[i]=Y1[i]+R.at（2,1）*（Z[i]-（*Zs））;
}
float*Z1=newfloat[total_number];
for（inti=0;i{
Z1[i]=0;
Z1[i]=Z1[i]+R.at（0,2）*（X[i]-（*Xs））;
Z1[i]=Z1[i]+R.at（1,2）*（Y[i]-（*Ys））;
Z1[i]=Z1[i]+R.at（2,2）*（Z[i]-（*Zs））;
}
//计算V=Ax-l中的A矩阵
cv:
:
MatA=Mat（2*total_number,6,CV_32F）;
for（inti=0;i<2*total_number;i++）
{
for（intj=0;j<6;j++）
{
A.at（i,j）=0;
}
}
//对线元素赋值
for（inti=0;i<2*total_number;i=i+2）
{
A.at（i+0,0）=（1.0/Z1[i/2]）*（R.at（0,0）*fk+R.at（0,2）*（x[i/2]-x00））;
A.at（i+1,0）=（1.0/Z1[i/2]）*（R.at（0,1）*fk+R.at（0,2）*（y[i/2]-y00））;
A.at（i+0,1）=（1.0/Z1[i/2]）*（R.at（1,0）*fk+R.at（1,2）*（x[i/2]-x00））;
A.at（i+1,1）=（1.0/Z1[i/2]）*（R.at（1,1）*fk+R.at（1,2）*（y[i/2]-y00））;
A.at（i+0,2）=（1.0/Z1[i/2]）*（R.at（2,0）*fk+R.at（2,2）*（x[i/2]-x00））;
A.at（i+1,2）=（1.0/Z1[i/2]）*（R.at（2,1）*fk+R.at（2,2）*（y[i/2]-y00））;
}
//对角元素赋值
for（inti=0;i<2*total_number;i=i+2）
{
A.at（i+0,3）=+1*（y[i/2]-y00）*sin（*w）-（（（x[i/2]-x00）/fk）*（（x[i/2]-x00）*cos（*k）-（y[i/2]-y00）*sin（*k））+fk*cos（*k））*cos（*w）;
A.at（i+1,3）=-1*（x[i/2]-x00）*sin（*w）-（（（y[i/2]-y00）/fk）*（（x[i/2]-x00）*cos（*k）-（y[i/2]-y00）*sin（*k））-fk*sin（*k））*cos（*w）;
A.at（i+0,4）=-1*fk*sin（*k）-（（x[i/2]-x00）/fk）*（（x[i/2]-x00）*sin（*k）+（y[i/2]-y00）*cos（*k））;
A.at（i+1,4）=-1*fk*cos（*k）-（（y[i/2]-y00）/fk）*（（x[i/2]-x00）*sin（*k）+（y[i/2]-y00）*cos（*k））;
A.at（i+0,5）=+1*（y[i/2]-y00）;
A.at（i+1,5）=-1*（x[i/2]-x00）;
}
//计算V=Ax-l中的l矩阵
cv:
:
Matl=Mat（2*total_number,1,CV_32F）;
for（inti=0;i<2*total_number;i=i+2）
{
l.at（i+0,0）=（x[i/2]-（-1*fk*X1[i/2]/Z1[i/2]））;
l.at（i+1,0）=（y[i/2]-（-1*fk*Y1[i/2]/Z1[i/2]））;
}
//已知Al矩阵之后，根据V=Ax-l来计算x
cv:
:
Matdx=Mat（6,1,CV_32F）
升级会员 