摄影测量作业3-空间后方交会计算

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);
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; }
4
10.46 64.43 40426.54 30319.81 757.31
以单像空间后方交会方法，求解该像片的外方位元素。
作业要求： 1、用一门计算机语言（如 C，C++，C#，VB）编写单像空间后方交会程序，各角
元素迭代计算至其改正值小于 6 秒。 2、提交正式的课程作业报告。 3、作业报告包括：封面、目录、正文等，其中正文部分包括：作业任务、计算
地球科学与环境工程学院 摄影测量实验报告书
课 程 名： 《摄影测量学》 学 号： 姓 名： 指导老师： 日 期： 2017 年 4 月 20 日
1
目录
一、目的与要求 ................................................ 3 二、实验数据与实验准备 ........................................ 3 三、程序流程图 ................................................ 4 四、实验代码 .................................................. 5 五、实验结果截图 ............................................. 10 六、实验总结 ................................................. 11
原理、算法流程、源程序、计算结果、结果分源自文库、心得体会等。
三．实验所用到的数学公式及程序计算步骤。
单张影像的空间后方交会：利用已知地面控制点数据及相应像点坐标 根据共线方程反 求影像的外方位元素。 数学模型：共线条件方程式：
3
求解过程： （1）获取已知数据。从航摄资料中查取平均航高与摄影机主距；获取控制点的地面测
int 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;
点号
像点坐标
地面坐标
x(mm) y(mm) X(m)
Y(m)
Z(m)
1
-86.15 -68.99 36589.41 25273.32 2195.17
2
-53.40 82.21 37631.08 31324.51 728.69
3
-14.78 -76.63 39100.97 24934.98 2386.50
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);
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)//计算系数矩 阵A {
2
一．实验目的
掌握摄影测量空间后方交会的原理，利用计算机编程语言实现空间后方交 会外方位元素的解算。
二．实验数据及实验设备
计算机 windows10 系统，编程软件 VS2015，地面控制点在摄影测量坐标系 中的坐标及其像点坐标文件 data.txt。
已知条件
摄影机主距 f=153.24mm，x0=0.01mm，y0=-0.02mm, 像片比例尺为 1:40000，像片上四个像点的框标坐标及其对应的地面坐标如下表。
m_list1.SetItemText(index1, i + 1, strSplit[i]);
} }
} sf.Close();
6
int Line = m_list1.GetItemCount(); CMatrix listXYZ(Line, 5); for (int i = 0; i < Line; i++)
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;
//以读的形式打开文件，如果打开失败 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) {
for (int i = 0; i < iRow; i++)
7
{ 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;
5
CFileDialog dlgOpenFile(TRUE, _T("txt"), NULL, OFN_FILEMUSTEXIST, _T("(文本文件)|*.txt|(所有文件）|*.*)||"));
if (dlgOpenFile.DoModal() == IDCANCEL) return;//如果选择取消按钮，则退出
量坐标并转换为地面摄影测量坐标。 （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），直到满足要求为止。
四．程序流程图。
4
五．程序的主要源代码如下所示：
#include "stdafx.h" #include "KongJianHouFangJiaoHui.h" #include "KongJianHouFangJiaoHuiDlg.h" #include "afxdialogex.h" #include "CeHuiLei.h" #include "Matrix.h" #ifdef _DEBUG #define new DEBUG_NEW #endif CMatrix xyXYZ; int ButtonClickeOk = 1; void CKongJianHouFangJiaoHuiDlg::OnBnClickedOk()//导入测量数据按钮代码 {
CString strPathName = dlgOpenFile.GetPathName();//获取选择的文件的完整路径
CString strFileTitle = dlgOpenFile.GetFileTitle();//获取文件名 CString strExtName = dlgOpenFile.GetFileExt();//获取文件扩展名 CStdioFile sf; //创建文件对象
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++) {
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)//计算外方位元素改正数矩 阵 {