摄影测量作业3.docx

地球科学与环境工程学院摄影测量实验报告书 课程名：摄影测量测量学 学 号： 2014113248 姓 名： 黄黎东 指导老师： 陈 强 日 期： 2016年4月17日18目录一、目的与要求3二、实验数学模型及算法分析3三、程序流程图4四、实验代码5五、实验结果截图10六、实验总结11一实验目的 掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交会外方位元素的解算。二实验数据及实验设备 计算机windows7系统，编程软件VS2015，地面控制点在摄影测量坐标系中的坐标及其像点坐标文件data.txt。已知条件摄影机主距f=153.24mm，x0=0.01mm，y0=-0.02mm, 像片比例尺为1:40000，像片上四个像点的框标坐标及其对应的地面坐标如下表。 点号像点坐标地面坐标x(mm)y(mm)X(m)Y(m)Z(m)1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31以单像空间后方交会方法，求解该像片的外方位元素。三实验所用到的数学公式及程序计算步骤。单张影像的空间后方交会：利用已知地面控制点数据及相应像点坐标 根据共线方程反求影像的外方位元素。 数学模型：共线条件方程式：求解过程：（1）获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测量坐标并转换为地面摄影测量坐标。（2）量测控制点的像点坐标并做系统改正。（3）确定未知数的初始值。在竖直摄影且地面控制点大致分布均匀的情况下，按如下方法确定初始值，即： q=w=k=0 式中:m为摄影比例尺分母；n为控制点个数。 （4）用三个角元素的初始值，计算个方向余弦，组成旋转矩阵R。（5）逐点计算像点坐标的近似值。利用未知数的近似值和控制点的地面坐标代入共线方程式，逐点计算像点坐标的近似值（x）、（y）。（6）逐点计算误差方程式的系数和常数项，组成误差方程式。（7）计算法方程的系数矩阵A和常数项L，组成法方程式。（8）解法方程，求得外方位元素的改正数dXs，dYs，dZs，dq，dw，dk。（9）用前次迭代取得的近似值，加本次迭代的改正数，计算外方位元素的新值。（10） 将求得的外方位元素改正数与规定的限差比较，若小于限差则迭代结束。否则用新的近似值重复（4）（9），直到满足要求为止。四程序流程图。 五程序的主要源代码如下所示： #include stdafx.h#include KongJianHouFangJiaoHui.h#include KongJianHouFangJiaoHuiDlg.h#include afxdialogex.h#include CeHuiLei.h#include Matrix.h#ifdef _DEBUG#define new DEBUG_NEW#endifCMatrix xyXYZ;int ButtonClickeOk = 1;void CKongJianHouFangJiaoHuiDlg:OnBnClickedOk()/导入测量数据按钮代码CFileDialog dlgOpenFile(TRUE, _T(txt), NULL,OFN_FILEMUSTEXIST, _T(文本文件)|*.txt|(所有文件）|*.*)|);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;CeHuiLei SplitString;int hang = 0;while (sf.ReadString(strLine)if (hang = 1)hang+;else int n = SplitString.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 5; i+)m_list1.SetItemText(index1, i + 1, strSpliti);sf.Close();int Line = m_list1.GetItemCount();CMatrix listXYZ(Line, 5);for (int i = 0; i Line; i+)for (int j = 0; j 5;j+)listXYZ(i,j)=_wtof(m_list1.GetItemText(i, j + 1);xyXYZ = listXYZ;ButtonClickeOk = 0;/表示导入数据按钮已被点击并且已导入数据CMatrix CKongJianHouFangJiaoHuiDlg:GetX(CMatrix A, CMatrix L)/计算外方位元素改正数矩阵CMatrix X,_A,_AA,N_AA;_A = A;/A的转置_AA = _A*A;N_AA = _AA.Inv();/_AA的逆矩阵X = N_AA*_A*L;return X;CMatrix CKongJianHouFangJiaoHuiDlg:GetA(CMatrix xyXYZ, double f, CMatrix XX)/计算系数矩阵Aint iRow = xyXYZ.Row();CMatrix A(2 * iRow, 6);double XS = XX(0, 0);double YS = XX(0, 1);double ZS = XX(0, 2);double w = XX(0, 4);double k = XX(0, 5);double q = XX(0, 3);double a1 = cos(q)*cos(k) - sin(q)*sin(w)*sin(k);double a2 = -cos(q)*sin(k) - sin(q)*sin(w)*cos(k);double a3 = -sin(q)*cos(w);double b1 = cos(w)*sin(k);double b2 = cos(w)*cos(k);double b3 = -sin(w);double c1 = sin(q)*cos(k) + cos(q)*sin(w)*sin(k);double c2 = -sin(q)*sin(k) + cos(q)*sin(w)*cos(k);double c3 = cos(q)*cos(w);f = f / 1000.0;for (int i = 0; i iRow; i+)double ZA = xyXYZ(i, 4);double XA = xyXYZ(i, 2);double YA = xyXYZ(i, 3);double x = (xyXYZ(i, 0)-_wtof(strx0)/1000.0;double y = (xyXYZ(i, 1)-_wtof(stry0)/1000.0;double _Z = a3*(XA - XS) + b3*(YA - YS) + c3*(ZA - ZS);A(2*i, 0) = (a1*f + a3*x) / _Z;A(2*i, 1) = (b1*f + b3*x) / _Z;A(2*i, 2) = (c1*f + c3*x) / _Z;A(2*i, 3) = y*sin(w) - (x*(x*cos(k) - y*sin(k) / f + f*cos(k)*cos(w);A(2*i, 4) = -f*sin(k) - x*(x*sin(k) + y*cos(k) / f;A(2*i, 5) = y;A(2*i+1, 0) = (a2*f + a3*y) / _Z;A(2 * i + 1, 1) = (b2*f + b3*y) / _Z;A(2 * i + 1, 2) = (c2*f + c3*y) / _Z;A(2 * i + 1, 3) = -x*sin(w) - (y*(x*cos(k) - y*sin(k) / f - f*sin(k)*cos(w);A(2 * i + 1, 4) = -f*cos(k) - y/ f*(x*sin(k) + y*cos(k);A(2 * i + 1, 5) = -x;return A;CMatrix CKongJianHouFangJiaoHuiDlg:GetL(CMatrix xyXYZ, double f, CMatrix XX)/计算L矩阵int iRow = xyXYZ.Row();CMatrix L(2 * iRow, 1);double XS = XX(0, 0);double YS = XX(0, 1);double ZS = XX(0, 2);double w = XX(0, 4);double k = XX(0, 5);double q = XX(0, 3);double a1 = cos(q)*cos(k) - sin(q)*sin(w)*sin(k);double a2 = -cos(q)*sin(k) - sin(q)*sin(w)*cos(k);double a3 = -sin(q)*cos(w);double b1 = cos(w)*sin(k);double b2 = cos(w)*cos(k);double b3 = -sin(w);double c1 = sin(q)*cos(k) + cos(q)*sin(w)*sin(k);double c2 = -sin(q)*sin(k) + cos(q)*sin(w)*cos(k);double c3 = cos(q)*cos(w);f = f / 1000.0;for (int i = 0; i iRow; i+)double ZA = xyXYZ(i, 4);double XA = xyXYZ(i, 2);double YA = xyXYZ(i, 3);double x = (xyXYZ(i, 0) - _wtof(strx0) / 1000.0;double y = (xyXYZ(i, 1) - _wtof(stry0) / 1000.0;double x0 = -f*(a1*(XA - XS) + b1*(YA - YS) + c1*(ZA - ZS) / (a3*(XA - XS) + b3*(YA - YS) + c3*(ZA - ZS);double y0 = -f*(a2*(XA - XS) + b2*(YA - YS) + c2*(ZA - ZS) / (a3*(XA - XS) + b3*(YA - YS) + c3*(ZA - ZS);L(2*i, 0) =x - x0;L(2*i+1, 0) =y - y0;return L;void CKongJianHouFangJiaoHuiDlg:OnBnClickedButton1()/计算按钮UpdateData(true);if (ButtonClickeOk = 1)MessageBox(_T(请先点击“导入已知坐标”按钮导入数据！), _T(提示);return;else if (strf = _T() | strm = _T() | strx0 = _T() | stry0 = _T()MessageBox(_T(请输入数据！), _T(提示);return;elseint c = strm.GetLength();double m = _wtof(strm.Right(c - 2);int iRow = xyXYZ.Row();double sumX = 0;double sumY = 0;double sumZ = 0;for (int i = 0; i iRow; i+)sumX += xyXYZ(i, 2);sumY += xyXYZ(i, 3);sumZ += xyXYZ(i, 4);double XS = sumX / iRow;double YS = sumY / iRow;double ZS = sumZ / iRow + m*_wtof(strf) / 1000;double w = 0;double k = 0;double q = 0;CMatrix XX(1, 6);XX(0, 0) = XS;XX(0, 1) = YS;XX(0, 2) = ZS;XX(0, 3) = q;XX(0, 4) = w;XX(0, 5) = k;CMatrix dX(6, 1);/迭代计算：do CMatrix A = GetA(xyXYZ, _wtof(strf), XX);CMatrix L = GetL(xyXYZ, _wtof(strf), XX);dX = GetX(A, L);for (int n = 0; n = 2.