摄影测量作业

1地球科学与环境工程学院摄 影 测 量 课 后 作 业 4课 程 名 ：摄影测量测量学学 号： 2014113248 姓 名： 黄黎东 指导老师： 陈 强 日 期： 2016 年 5 月 21 日2目录一、目的与要求 .3二、实验数学模型及算法分析 .3三、程序流程图 .5四、实验代码 .5五、实验结果截图 .13六、实验总结 .143一 实验目的 （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，计算公式如下：4计算两张影像各自测量得的框标坐标所对应的地面坐标，计算公式为：（2）共线方程严密法计算过程：根据共线条件方程推导基于共线严密方程的前方交会公式；共线条件方程为：根据共线条件方程推出共线条件严密方程形式如下：L1 X + L2 Y +L3 Z - Lx = 0L4 X + L5 Y +L6 Z - Ly = 0求出共线严密方程的系数 L1 ,L2, L3, L4, L5, L6, Lx, Ly.L1=fa1+(x-x0)a3 L2= fb1+(x-x0)b3 L3= fc1+(x-x0)c3L4= fa2+(y-y0)a3 L5= fb2+(y-y0)b3 L6= fc2+(y-y0)c3Lx=fa1Xs+fb1Ys+fc1Zs+(x-x0)a3Xs+(x-x0)b3Ys+(x-x0)c3ZsLy=fa2Xs+fb2Ys+fc2Zs+(y-y0)a3Xs+(y-y0)b3Ys+(y-y0)c3Zs用平差方法计算出地面坐标 根据严密方程得到系数阵为：A= L=123456 根据最小二乘间接平差原理推导出如下式子：(ATPA) =ATP 则 =(ATA)-1ATL四程序流程图。5五程序的主要源代码如下所示：#include “stdafx.h“#include “HomeWork4.h“#include “HomeWork4Dlg.h“#include “afxdialogex.h“#include “Matrix.h“/矩阵类#ifdef _DEBUG#define new DEBUG_NEW#endifCMatrix xy1, xy2, XYZqwk(2,6);bool ButtonOpenFile = false;BOOL CHomeWork4Dlg:OnInitDialog()CDialogEx:OnInitDialog();strx0 = _T(“-0.004“);stry0 = _T(“-0.008“);strf = _T(“152.72“);开始读取文件点投影法 or严密法？点投影法严密方程法点投影法数学模型计算处理严密方程法数学模型计算处理输出计算结果结束6UpdateData(false);DWORD dwStyle = 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);DWORD dwStyle1 = 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 ASSERT(IDM_ABOUTBOX AppendMenu(MF_SEPARATOR);pSysMenu-AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);/ 设置此对话框的图标。 当应用程序主窗口不是对话框时，框架将自动/ 执行此操作7SetIcon(m_hIcon, TRUE); / 设置大图标SetIcon(m_hIcon, FALSE); / 设置小图标/ TODO: 在此添加额外的初始化代码return TRUE; / 除非将焦点设置到控件，否则返回 TRUEvoid CHomeWork4Dlg:OnBnClickedOk()/读取数据按钮代码CFileDialog dlgOpenFile(TRUE, _T(“csv“), NULL,OFN_FILEMUSTEXIST, _T(“(文本文件)|*.csv|(所有文件）|*.*)|“);if (dlgOpenFile.DoModal() = IDCANCEL)return;/如果选择取消按钮，则退出CString strPathName = dlgOpenFile.GetPathName();/获取选择的文件的完整路径CString strFileTitle = dlgOpenFile.GetFileTitle();/获取文件名CString strExtName = dlgOpenFile.GetFileExt();/获取文件扩展名CStdioFile sf; /创建文件对象/以读的形式打开文件，如果打开失败if (!sf.Open(strPathName, CFile:modeRead)MessageBox(_T(“读取文件出错！“);return;CString strLine;/存放每一行文本CStringArray strSplit;int hang = 0;while (sf.ReadString(strLine)if (hang =0)hang+;else if(hang=1|hang=2)int n =SplitStringArray(strLine, , strSplit);int index1 = m_list1.GetItemCount();CString sno1 = _T(“);sno1.Format(_T(“%d“), index1 + 1);int nR = m_list1.InsertItem(index1, sno1);for (int i = 0; i 2)m_list1.SetItemText(index1, i + 2, strSpliti);elsem_list1.SetItemText(index1, i + 1, strSpliti);int index2 = m_list2.GetItemCount();CString sno2 = _T(“);if(hang=1)sno2=_T(“左片“);else sno2= _T(“右片“);int nR2 = m_list2.InsertItem(index2, sno2);for (int j = 4; j 2)m_list1.SetItemText(index1, i + 2, strSpliti);elsem_list1.SetItemText(index1, i + 1, strSpliti);9sf.Close();ButtonOpenFile =true;int CHomeWork4Dlg:SplitStringArray(CString str, char split, CStringArrayint idx = str.Find(split, startIdx);aStr.RemoveAll();/先清空while (-1 != idx)CString sTmp = str.Mid(startIdx, idx - startIdx);aStr.Add(sTmp);startIdx = idx + 1;idx = str.Find(split, startIdx);CString sTmp = str.Right(str.GetLength() - startIdx);if (!sTmp.IsEmpty()aStr.Add(sTmp);return aStr.GetSize();CMatrix CHomeWork4Dlg:XYZ(CMatrix XYZqwk, CMatrix xy1, CMatrix xy2,double f)/点投影系数法函数模型double PI = 3.141592654;double Bx = XYZqwk(1, 0) - XYZqwk(0, 0);double Bz = XYZqwk(1, 2) - XYZqwk(0, 2);CMatrix R1(3, 3),R2(3,3);double q = XYZqwk(0, 3)*PI / 180;double w= XYZqwk(0, 4)*PI / 180;double k = XYZqwk(0, 5)*PI / 180;double q1 = XYZqwk(1, 3)*PI / 180;double w1 = XYZqwk(1, 4)*PI / 180;double k1 = 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);10R2(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;int iRow = xy1.Row();CMatrix xyf1(3, iRow),xyf2(3,iRow),XYZ1(3,iRow),XYZ2(3,iRow),XYZ(3,iRow);for (int i = 0; i iRow; i+)xyf1(0, i) = xy1(i, 0)/1000;