摄影测量作业.docx
《摄影测量作业.docx》由会员分享，可在线阅读，更多相关《摄影测量作业.docx（16页珍藏版）》请在冰豆网上搜索。
摄影测量作业
地球科学与环境工程学院
摄影测量课后作业4
课程名：
《摄影测量测量学》
学号：
2014113248
姓名：
黄黎东
指导老师：
陈强
日期：
2016年5月21日
一、目的与要求3
二、实验数学模型及算法分析3
三、程序流程图5
四、实验代码5
五、实验结果截图13
六、实验总结14
一．实验目的
（1）掌握摄影测量立体像对的空间前方交会的计算方法及过程；
（2）比较采用点投影系数法和共线方程严密法计算地面点坐标的差异。
二．实验数据及实验设备
计算机windows10系统，编程软件VS2015，两张影像各自的外方位元素和框标坐标及相机的内方位元素存在文件data.excel中。
已知条件
左像片的外方位元素为:
（14922.9778,11532.4077,3230.3383）m,（-0.0566,-0.1215,-0.6456）°，
右像片的外方位元素为:
（16296.9601,11524.2256,3239.6070）m,（0.0429,-0.5148,-0.1538）°，
相机的内方位元素（x0,y0,f）为（-0.004，-0.008,152.72）mm。
现量测了6个像点的框标坐标，如下（单位mm）：
左片:
1（85.7255,69.6561）,2（0.3206,59.1841）,3（83.7954,-1.9617）,
4（-2.2603,-15.1519）,5（88.2371,-73.7021）,6（2.4792,-67.5320）
右片:
1（0.6313,70.1653）,2（-84.4506,60.4553）,3（-3.3174,-1.3205）
4（-97.0211,-13.7425）,5（2.7516,-72.5590）,6（-82.9161,-65.7818）
采用点投影系数法和共线方程严密法，分别计算上述6个点的地面坐标。
三．实验所用到的数学公式及程序计算步骤。
（1）点投影系数法计算过程：
根据两张影像量测得的框标坐标用旋转矩阵分别计算其对应的像空间辅助坐标，计算公式如下：
根据两张影像的外方位元素计算摄影基线的三个分量，计算公式如下：
计算两张影像的投影系数N1,N2，计算公式如下：
计算两张影像各自测量得的框标坐标所对应的地面坐标，计算公式为：
（2）共线方程严密法计算过程：
根据共线条件方程推导基于共线严密方程的前方交会公式；
共线条件方程为：
根据共线条件方程推出共线条件严密方程形式如下：
L1X+L2Y+L3Z-Lx=0
L4X+L5Y+L6Z-Ly=0
求出共线严密方程的系数L1,L2,L3,L4,L5,L6,Lx,Ly.
L1=fa1+（x-x0）a3L2=fb1+（x-x0）b3L3=fc1+（x-x0）c3
L4=fa2+（y-y0）a3L5=fb2+（y-y0）b3L6=fc2+（y-y0）c3
Lx=fa1Xs+fb1Ys+fc1Zs+（x-x0）a3Xs+（x-x0）b3Ys+（x-x0）c3Zs
Ly=fa2Xs+fb2Ys+fc2Zs+（y-y0）a3Xs+（y-y0）b3Ys+（y-y0）c3Zs
用平差方法计算出地面坐标
根据严密方程得到系数阵为：
A=
L=
根据最小二乘间接平差原理推导出如下式子：
（ATPA）
=ATP
则
=（ATA）-1ATL
四．程序流程图。
五．程序的主要源代码如下所示：
#include"stdafx.h"
#include"HomeWork4.h"
#include"HomeWork4Dlg.h"
#include"afxdialogex.h"
#include"Matrix.h"//矩阵类
#ifdef_DEBUG
#definenewDEBUG_NEW
#endif
CMatrixxy1,xy2,XYZqwk（2,6）;
boolButtonOpenFile=false;
BOOLCHomeWork4Dlg:
:
OnInitDialog（）
{
CDialogEx:
:
OnInitDialog（）;
strx0=_T（"-0.004"）;
stry0=_T（"-0.008"）;
strf=_T（"152.72"）;
UpdateData（false）;
DWORDdwStyle=m_list1.GetExtendedStyle（）;
dwStyle|=LVS_EX_FULLROWSELECT;//选中某行使整行高亮（只适用与report风格的listctrl）
dwStyle|=LVS_EX_GRIDLINES;//网格线（只适用与report风格的listctrl）
m_list1.SetExtendedStyle（dwStyle）;//设置扩展风格
m_list1.InsertColumn（0,_T（"左片"）,LVCFMT_CENTER,40）;
m_list1.InsertColumn（1,_T（"x（mm）"）,LVCFMT_CENTER,100）;
m_list1.InsertColumn（2,_T（"y（mm）"）,LVCFMT_CENTER,100）;
m_list1.InsertColumn（3,_T（"右片"）,LVCFMT_CENTER,40）;
m_list1.InsertColumn（4,_T（"x（mm）"）,LVCFMT_CENTER,100）;
m_list1.InsertColumn（5,_T（"y（mm）"）,LVCFMT_CENTER,100）;
DWORDdwStyle1=m_list2.GetExtendedStyle（）;
dwStyle1|=LVS_EX_FULLROWSELECT;//选中某行使整行高亮（只适用与report风格的listctrl）
dwStyle1|=LVS_EX_GRIDLINES;//网格线（只适用与report风格的listctrl）
m_list2.SetExtendedStyle（dwStyle1）;//设置扩展风格
m_list2.InsertColumn（0,_T（""）,LVCFMT_CENTER,40）;
m_list2.InsertColumn（1,_T（"Xs（m）"）,LVCFMT_CENTER,100）;
m_list2.InsertColumn（2,_T（"Ys（m）"）,LVCFMT_CENTER,100）;
m_list2.InsertColumn（3,_T（"Zs（m）"）,LVCFMT_CENTER,100）;
m_list2.InsertColumn（4,_T（"Fi（度）"）,LVCFMT_CENTER,100）;
m_list2.InsertColumn（5,_T（"Omega（度）"）,LVCFMT_CENTER,100）;
m_list2.InsertColumn（6,_T（"Kama（度）"）,LVCFMT_CENTER,100）;
ASSERT（（IDM_ABOUTBOX&0xFFF0）==IDM_ABOUTBOX）;
ASSERT（IDM_ABOUTBOX<0xF000）;
CMenu*pSysMenu=GetSystemMenu（FALSE）;
if（pSysMenu!
=NULL）
{
BOOLbNameValid;
CStringstrAboutMenu;
bNameValid=strAboutMenu.LoadString（IDS_ABOUTBOX）;
ASSERT（bNameValid）;
if（!
strAboutMenu.IsEmpty（））
{
pSysMenu->AppendMenu（MF_SEPARATOR）;
pSysMenu->AppendMenu（MF_STRING,IDM_ABOUTBOX,strAboutMenu）;
}
}
//设置此对话框的图标。
当应用程序主窗口不是对话框时，框架将自动
//执行此操作
SetIcon（m_hIcon,TRUE）;//设置大图标
SetIcon（m_hIcon,FALSE）;//设置小图标
//TODO:
在此添加额外的初始化代码
returnTRUE;//除非将焦点设置到控件，否则返回TRUE}
voidCHomeWork4Dlg:
:
OnBnClickedOk（）//读取数据按钮代码
{
CFileDialogdlgOpenFile（TRUE,_T（"csv"）,NULL,
OFN_FILEMUSTEXIST,_T（"（文本文件）|*.csv|（所有文件）|*.*）||"））;
if（dlgOpenFile.DoModal（）==IDCANCEL）
return;//如果选择取消按钮，则退出
CStringstrPathName=dlgOpenFile.GetPathName（）;//获取选择的文件的完整路径
CStringstrFileTitle=dlgOpenFile.GetFileTitle（）;//获取文件名
CStringstrExtName=dlgOpenFile.GetFileExt（）;//获取文件扩展名
CStdioFilesf;//创建文件对象
//以读的形式打开文件，如果打开失败
if（!
sf.Open（strPathName,CFile:
:
modeRead））
{
MessageBox（_T（"读取文件出错！
"））;
return;
}
CStringstrLine;//存放每一行文本
CStringArraystrSplit;
inthang=0;
while（sf.ReadString（strLine））
{
if（hang==0）
{
hang++;
}
elseif（hang==1||hang==2）
{
intn=SplitStringArray（strLine,',',strSplit）;
intindex1=m_list1.GetItemCount（）;
CStringsno1=_T（""）;
sno1.Format（_T（"%d"）,index1+1）;
intnR=m_list1.InsertItem（index1,sno1）;
for（inti=0;i<4;i++）
{
if（i==2）
{
m_list1.SetItemText（index1,i+1,sno1）;
m_list1.SetItemText（index1,i+2,strSplit[i]）;
}
elseif（i>2）
{
m_list1.SetItemText（index1,i+2,strSplit[i]）;
}
else
m_list1.SetItemText（index1,i+1,strSplit[i]）;
}
intindex2=m_list2.GetItemCount（）;
CStringsno2=_T（""）;
if（hang==1）
sno2=_T（"左片"）;
else
sno2=_T（"右片"）;
intnR2=m_list2.InsertItem（index2,sno2）;
for（intj=4;j<10;j++）
{
m_list2.SetItemText（index2,j-3,strSplit[j]）;
}
hang++;
}
else
{
intn=SplitStringArray（strLine,',',strSplit）;
intindex1=m_list1.GetItemCount（）;
CStringsno1=_T（""）;
sno1.Format（_T（"%d"）,index1+1）;
intnR=m_list1.InsertItem（index1,sno1）;
for（inti=0;i<4;i++）
{
if（i==2）
{
m_list1.SetItemText（index1,i+1,sno1）;
m_list1.SetItemText（index1,i+2,strSplit[i]）;
}
elseif（i>2）
{
m_list1.SetItemText（index1,i+2,strSplit[i]）;
}
else
m_list1.SetItemText（index1,i+1,strSplit[i]）;
}
}
}
sf.Close（）;
ButtonOpenFile=true;
}
intCHomeWork4Dlg:
:
SplitStringArray（CStringstr,charsplit,CStringArray&aStr）//字符串切割成员函数代码
{
intstartIdx=0;
intidx=str.Find（split,startIdx）;
aStr.RemoveAll（）;//先清空
while（-1!
=idx）
{
CStringsTmp=str.Mid（startIdx,idx-startIdx）;
aStr.Add（sTmp）;
startIdx=idx+1;
idx=str.Find（split,startIdx）;
}
CStringsTmp=str.Right（str.GetLength（）-startIdx）;
if（!
sTmp.IsEmpty（））
aStr.Add（sTmp）;
returnaStr.GetSize（）;
}
CMatrixCHomeWork4Dlg:
:
XYZ（CMatrixXYZqwk,CMatrixxy1,CMatrixxy2,doublef）//点投影系数法函数模型
{
doublePI=3.141592654;
doubleBx=XYZqwk（1,0）-XYZqwk（0,0）;
doubleBz=XYZqwk（1,2）-XYZqwk（0,2）;
CMatrixR1（3,3）,R2（3,3）;
doubleq=XYZqwk（0,3）*PI/180;
doublew=XYZqwk（0,4）*PI/180;
doublek=XYZqwk（0,5）*PI/180;
doubleq1=XYZqwk（1,3）*PI/180;
doublew1=XYZqwk（1,4）*PI/180;
doublek1=XYZqwk（1,5）*PI/180;
R1（0,0）=cos（q）*cos（k）-sin（q）*sin（w）*sin（k）;
R1（0,1）=-cos（q）*sin（k）-sin（q）*sin（w）*cos（k）;
R1（0,2）=-sin（q）*cos（w）;
R1（1,0）=cos（w）*sin（k）;
R1（1,1）=cos（w）*cos（k）;
R1（1,2）=-sin（w）;
R1（2,0）=sin（q）*cos（k）+cos（q）*sin（w）*sin（k）;
R1（2,1）=-sin（q）*sin（k）+cos（q）*sin（w）*cos（k）;
R1（2,2）=cos（q）*cos（w）;
R2（0,0）=cos（q1）*cos（k1）-sin（q1）*sin（w1）*sin（k1）;
R2（0,1）=-cos（q1）*sin（k1）-sin（q1）*sin（w1）*cos（k1）;
R2（0,2）=-sin（q1）*cos（w1）;
R2（1,0）=cos（w1）*sin（k1）;
R2（1,1）=cos（w1）*cos（k1）;
R2（1,2）=-sin（w1）;
R2（2,0）=sin（q1）*cos（k1）+cos（q1）*sin（w1）*sin（k1）;
R2（2,1）=-sin（q1）*sin（k1）+cos（q1）*sin（w1）*cos（k1）;
R2（2,2）=cos（q1）*cos（w1）;
f=f/1000;
intiRow=xy1.Row（）;
CMatrixxyf1（3,iRow）,xyf2（3,iRow）,XYZ1（3,iRow）,XYZ2（3,iRow）,XYZ（3,iRow）;
for（inti=0;i{
xyf1（0,i）=xy1（i,0）/1000;
xyf1（1,i）=xy1（i,1）/1000;
xyf1（2,i）=-f;
xyf2（0,i）=xy2（i,0）/1000;
xyf2（1,i）=xy2（i,1）/1000;
xyf2（2,i）=-f;
}
XYZ1=R1*xyf1;
XYZ2=R2*xyf2;
doubleN1;
for（intj=0;j{
N1=（Bx*XYZ2（2,j）-Bz*XYZ2（0,j））/（XYZ1（0,j）*XYZ2（2,j）-XYZ2（0,j）*XYZ1（2,j））;
XYZ（0,j）=XYZqwk（0,0）+N1*XYZ1（0,j）;
XYZ（1,j）=XYZqwk（0,1）+N1*XYZ1（1,j）;
XYZ（2,j）=XYZqwk（0,2）+N1*XYZ1（2,j）;
}
returnXYZ;
}
CMatrixCHomeWork4Dlg:
:
inlineXYZ（CMatrixXYZqwk,CMatrixxy1,CMatrixxy2,doublef）//共线方程严密法函数
{
intn=xy1.Row（）;
CMatrixA（4,3）,L（4,1）,X,XYZ（3,n）;
doublePI=3.141592654;
doubleq[2],w[2],k[2];
doublea1[2],a2[2],a3[2],b1[2],b2[2],b3[2],c1[2],c2[2],c3[2];
for（inti=0;i<2;i++）
{
q[i]=XYZqwk（i,3）*PI/180;
w[i]=XYZqwk（i,4）*PI/180;
k[i]=XYZqwk（i,5）*PI/180;
a1[i]=cos（q[i]）*cos（k[i]）-sin（q[i]）*sin（w[i]）*sin（k[i]）;
a2[i]=-cos（q[i]）*sin（k[i]）-sin（q[i]）*sin（w[i]）*cos（k[i]）;
a3[i]=-sin（q[i]）*cos（w[i]）;
b1[i]=cos（w[i]）*sin（k[i]）;
b2[i]=cos（w[i]）*cos（k[i]）;
b3[i]=-sin（w[i]）;
c1[i]=sin（q[i]）*cos（k[i]）+cos（q[i]）*sin（w[i]）*sin（k[i]）;
c2[i]=-sin（q[i]）*sin（k[i]）+cos（q[i]）*sin（w[i]）*cos（k[i]）;
c3[i]=cos（q[i]）*cos（w[i]）;
}
f=f/1000.0;
for（intj=0;j{
A（0,0）=f*a1[0]+xy1（j,0）/1000.0*a3[0];
A（0,1）=f*b1[0]+xy1（j,0）/1000.0*b3[0];
A（0,2）=f*c1[0]+xy1（j,0）/1000.0*c3[0];
A（1,0）=f*a2[0]+xy1（j,1）/1000.0*a3[0];
A（1,1）=f*b2[0]+xy1（j,1）/1000.0*b3[0];
A（1,2）=f*c2[0]+xy1（j,1）/1000.0*c3[0];
A（2,0）=f*a1[1]+xy2（j,0）/1000.0*a3[1];
A（2,1）=f*b1[1]+xy2（j,0）/1000.0*b3[1];
A（2,2）=f*c1[1]+xy2（j,0）/1000.0*c3[1];
A（3,0）=f*a2[1]+xy2（j,1）/1000.0*a3[1];
A（3,1）=f*b2[1]+xy2（j,1）/1000.0*b3[1];
A（3,2）=f*c2[1]+xy2（j,1）/1000.0*c3[1];
L（0,0）=f*a1[0]*XYZqwk（0,0）+f*b1[0]*XYZqwk（0,1）+f*c1[0]*XYZqwk（0,2）+xy1（j,0）/1000.0*a3[0]*XYZqwk（0,0）+xy1（j,0）/1000.0*b3[0]*XYZqwk（0,1）+xy1（j,0）/1000.0*c3[0]*XYZqwk（0,2）;
L（1,0）=f*a2[0]*XYZqwk（0,0
升级会员 