《摄影测量测量学》摄影测量作业4.docx

地球科学与环境工程学院摄影测量课后作业4课程名：摄影测量测量学学 号： 姓 名： 指导老师： 日 期： 年 月 日目录一、目的与要求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）共线方程严密法计算过程：根据共线条件方程推导基于共线严密方程的前方交会公式； 共线条件方程为：根据共线条件方程推出共线条件严密方程形式如下：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=L1L2L3L4L5L6 L=LxLy根据最小二乘间接平差原理推导出如下式子：(ATPA)XYZ=ATPL则 XYZ=(ATA)-1ATL四程序流程图。开始读取文件点投影法or严密法？严密方程法点投影法点投影法数学模型计算处理严密方程法数学模型计算处理输出计算结果结束五程序的主要源代码如下所示：#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);UpdateData(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 & 0xFFF0) = IDM_ABOUTBOX);ASSERT(IDM_ABOUTBOX AppendMenu(MF_SEPARATOR);pSysMenu-AppendMenu(MF_STRING, IDM_ABOUTBOX, strAboutMenu);/ 设置此对话框的图标。 当应用程序主窗口不是对话框时，框架将自动/ 执行此操作SetIcon(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 10; j+)m_list2.SetItemText(index2, j-3, strSplitj);hang+;elseint 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);sf.Close();ButtonOpenFile =true;int CHomeWork4Dlg:SplitStringArray(CString str, char split, CStringArray& aStr)/字符串切割成员函数代码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()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);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;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;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+)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)/共线方程严密法函数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(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 = 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.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) + f*b20 * XYZqwk(0, 1) + f*c20 * XYZqwk(0, 2) + xy1(j, 1) / 1000.0 *a30 * XYZqwk(0, 0) + xy1(j, 1) / 1000.0 *b30 * XYZqwk(0, 1) + xy1(j, 1) / 1000.0 *c30 * XYZqwk(0, 2);L(2, 0) = f*a11 * XYZqwk(1, 0) + f*b11 * XYZqwk(1, 1) + f*c11 * XYZqwk(1, 2) + xy2(j, 0) / 1000.0 *a31 * XYZqwk(1, 0) + xy2(j, 0) / 1000.0 *b31 * XYZqwk(1, 1) + xy2(j, 0) / 1000.0 *c31 * XYZqwk(1, 2);L(3, 0) = f*a21 * XYZqwk(1, 0) + f*b21 * XYZqwk(1, 1) + f*c21 * XYZqwk(1, 2) + xy2(j, 1) / 1000.0 *a31 * XYZqwk(1, 0) + xy2(j, 1) / 1000.0 *b31 * XYZqwk(1, 1) + xy2(j, 1) / 1000.0 *c31 * XYZqwk(1, 2);CMatrix _A, ATA, N_AA;_A = A; ATA = _A*A;N_AA = ATA.Inv(); X= N_AA*_A*L; XYZ(0, j) = X(0, 0); XYZ(1, j) = X(1, 0); XYZ(2, j) = X(2, 0);return XYZ;void CHomeWork4Dlg:OnBnClickedButton1()/点投影系数法UpdateData(true);if (ButtonOpenFile = false)MessageBox(_T(请先导入数据！), _T(提示);return;else int Line = m_list1.GetItemCount();CMatrix listxy1(Line, 2), listxy2(Line, 2);for (int i = 0; i Line; i+)for (int j = 0; j 2; j+)if (j = 0)listxy1(i, j) = _wtof(m_list1.GetItemText(i, j + 1) - _wtof(strx0);listxy2(i, j) = _wtof(m_list1.GetItemText(i, j + 4) - _wtof(strx0);elselistxy1(i, j) = _wtof(m_list1.GetItemText(i, j + 1) - _wtof(stry0);listxy2(i, j) = _wtof(m_list1.GetItemText(i, j + 4) - _wtof(stry0);xy1 = listxy1;xy2 = listxy2;for (int k =