摄影测量作业.docx

1、摄影测量作业地球科学与环境工程学院摄影测量课后作业4课程名：摄影测量测量学学 号： 2014113248 姓 名： 黄黎东 指导老师： 陈 强 日 期： 2016年5月21日一、目的与要求 3二、实验数学模型及算法分析 3三、程序流程图 5四、实验代码 5五、实验结果截图 13六、实验总结 14一 实验目的 （1） 掌握摄影测量立体像对的空间前方交会的计算方法及过程；（2） 比较采用点投影系数法和共线方程严密法计算地面点坐标的差异。二实验数据及实验设备计算机windows10系统，编程软件VS2015，两张影像各自的外方位元素和框标坐标及相机的内方位元素存在文件data.excel中。已知条件

2、左像片的外方位元素为: (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,

3、 -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）点投影系数法计算过程：根据两张影像量测得的框标坐标用旋转矩阵分别计算其对应的像空间辅助坐标，计算公式如下：根据两张

4、影像的外方位元素计算摄影基线的三个分量，计算公式如下：计算两张影像的投影系数N1,N2，计算公式如下：计算两张影像各自测量得的框标坐标所对应的地面坐标，计算公式为：（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

5、= 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=根据最小二乘间接平差原理推导出如下式子：(ATPA)=ATP则 =(ATA)-1ATL四程序流程图。五程序的主要源代码如下所示：#include stdafx.h#include HomeWork4.h#include HomeW

6、ork4Dlg.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); UpdateData(false); DWORD dwStyl

7、e = 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, 10

8、0); 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;/

9、选中某行使整行高亮（只适用与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_lis

10、t2.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 AppendMenu(MF

11、_SEPARATOR); pSysMenu-AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu); / 设置此对话框的图标。 当应用程序主窗口不是对话框时，框架将自动 / 执行此操作 SetIcon(m_hIcon, TRUE); / 设置大图标 SetIcon(m_hIcon, FALSE); / 设置小图标 / TODO: 在此添加额外的初始化代码 return TRUE; / 除非将焦点设置到控件，否则返回 TRUEvoid CHomeWork4Dlg:OnBnClickedOk()/读取数据按钮代码 CFileDialog dlgOpenF

12、ile(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();/获

13、取文件扩展名 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

14、); 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); else m_list1.SetItemText(index1, i + 1, strSpliti); int index2 = m_list2.GetItemCount(); CStrin

15、g sno2 = _T(); if(hang=1) sno2=_T(左片); else sno2= _T(右片); int nR2 = m_list2.InsertItem(index2, sno2); for (int j = 4; j 10; j+) m_list2.SetItemText(index2, j-3, strSplitj); hang+; else int n = SplitStringArray(strLine, , strSplit); int index1 = m_list1.GetItemCount(); CString sno1 = _T(); sno1.Forma

16、t(_T(%d), index1 + 1); int nR = m_list1.InsertItem(index1, sno1); for (int i = 0; i 2) m_list1.SetItemText(index1, i + 2, strSpliti); else m_list1.SetItemText(index1, i + 1, strSpliti); sf.Close(); ButtonOpenFile =true;int CHomeWork4Dlg:SplitStringArray(CString str, char split, CStringArray& aStr)/字

17、符串切割成员函数代码 int startIdx = 0; int 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() a

18、Str.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

19、= 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) =

20、 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

21、, 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),xyf

22、2(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; 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; double N1; for (int j = 0; j iRow; j+

23、) 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); return XYZ;CMatrix CHomeWork4Dlg:inlineXYZ(CMatrix XYZqwk, CMatrix xy1, CMatrix xy2, double f)/共线方

24、程严密法函数 int n = xy1.Row(); CMatrix A(4,3),L(4,1),X,XYZ(3,n); double PI = 3.141592654; double q2, w2, k2; double a12, a22, a32, b12, b22, b32, c12, c22, c32; for (int i = 0; i 2; i+) qi = XYZqwk(i, 3)*PI / 180; wi = XYZqwk(i, 4)*PI / 180; ki = XYZqwk(i, 5)*PI / 180; a1i = cos(qi)*cos(ki) - sin(qi)*sin

25、(wi)*sin(ki); a2i = -cos(qi)*sin(ki) - sin(qi)*sin(wi)*cos(ki); a3i = -sin(qi)*cos(wi); b1i = cos(wi)*sin(ki); b2i = cos(wi)*cos(ki); b3i = -sin(wi); c1i = sin(qi)*cos(ki) + cos(qi)*sin(wi)*sin(ki); c2i = -sin(qi)*sin(ki) + cos(qi)*sin(wi)*cos(ki); c3i = cos(qi)*cos(wi); f = f / 1000.0; for (int j =

26、 0; j xy1.Row(); j+) A(0, 0) = f*a10 + xy1(j, 0)/1000.0*a30; A(0, 1) = f*b10 + xy1(j, 0)/1000.0 *b30; A(0, 2) = f*c10 + xy1(j, 0)/1000.0 *c30; A(1, 0) = f*a20 + xy1(j, 1)/1000.0 *a30; A(1, 1) = f*b20 + xy1(j, 1)/1000.0 *b30; A(1, 2) = f*c20 + xy1(j, 1)/1000.0 *c30; A(2, 0) = f*a11 + xy2(j, 0) / 1000

27、.0 *a31; A(2, 1) = f*b11 + xy2(j, 0) / 1000.0 *b31; A(2, 2) = f*c11 + xy2(j, 0) / 1000.0 *c31; A(3, 0) = f*a21 + xy2(j, 1) / 1000.0 *a31; A(3, 1) = f*b21 + xy2(j, 1) / 1000.0 *b31; A(3, 2) = f*c21 + xy2(j, 1) / 1000.0 *c31; L(0, 0) = f*a10 * XYZqwk(0, 0) + f*b10 * XYZqwk(0, 1) + f*c10 * XYZqwk(0, 2) + xy1(j, 0) / 1000.0 *a30 * XYZqwk(0, 0) + xy1(j, 0) / 1000.0 *b30 * XYZqwk(0, 1) + xy1(j, 0) / 1000.0 *c30 * XYZqwk(0, 2); L(1,0)= f*a20 * XYZqwk(0, 0