攝影測量作業(yè)3

1、地球科學(xué)與環(huán)境工程學(xué)院攝影測量實(shí)驗(yàn)報(bào)告書 課程名：攝影測量測量學(xué) 學(xué) 號(hào)： 2014113248 姓 名： 黃黎東 指導(dǎo)老師： 陳 強(qiáng) 日 期： 2016年4月17日目錄一、目的與要求3二、實(shí)驗(yàn)數(shù)學(xué)模型及算法分析3三、程序流程圖4四、實(shí)驗(yàn)代碼5五、實(shí)驗(yàn)結(jié)果截圖10六、實(shí)驗(yàn)總結(jié)11一實(shí)驗(yàn)?zāi)康?掌握攝影測量空間后方交會(huì)的原理，利用計(jì)算機(jī)編程語言實(shí)現(xiàn)空間后方交會(huì)外方位元素的解算。二實(shí)驗(yàn)數(shù)據(jù)及實(shí)驗(yàn)設(shè)備 計(jì)算機(jī)windows7系統(tǒng)，編程軟件VS2015，地面控制點(diǎn)在攝影測量坐標(biāo)系中的坐標(biāo)及其像點(diǎn)坐標(biāo)文件data.txt。已知條件攝影機(jī)主距f=153.24mm，x0=0.01mm，y0=-0.02mm,

2、像片比例尺為1:40000，像片上四個(gè)像點(diǎn)的框標(biāo)坐標(biāo)及其對(duì)應(yīng)的地面坐標(biāo)如下表。 點(diǎn)號(hào)像點(diǎn)坐標(biāo)地面坐標(biāo)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以單像空間后方交會(huì)方法，求解該像片的外方位元素。三實(shí)驗(yàn)所用到的數(shù)學(xué)公式及程序計(jì)算步驟。單張影像的空間后方交會(huì)：利用已知地面控制點(diǎn)數(shù)據(jù)及相應(yīng)像點(diǎn)坐標(biāo) 根據(jù)共線方程反求影像的外方位元

3、素。 數(shù)學(xué)模型：共線條件方程式：求解過程： （1）獲取已知數(shù)據(jù)。從航攝資料中查取平均航高與攝影機(jī)主距；獲取控制點(diǎn)的地面測量坐標(biāo)并轉(zhuǎn)換為地面攝影測量坐標(biāo)。 （2）量測控制點(diǎn)的像點(diǎn)坐標(biāo)并做系統(tǒng)改正。 （3）確定未知數(shù)的初始值。在豎直攝影且地面控制點(diǎn)大致分布均勻的情況下，按如下方法確定初始值，即： q=w=k=0  式中:m為攝影比例尺分母；n為控制點(diǎn)個(gè)數(shù)。 （4）用三個(gè)角元素的初始值，計(jì)算個(gè)方向余弦，組成旋轉(zhuǎn)矩陣R。 （5）逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值。利用未知數(shù)的近似值和控制點(diǎn)的地面坐標(biāo)代入共線方程式，逐點(diǎn)計(jì)算像點(diǎn)坐標(biāo)的近似值（x）、（y）。

4、0;（6）逐點(diǎn)計(jì)算誤差方程式的系數(shù)和常數(shù)項(xiàng)，組成誤差方程式。 （7）計(jì)算法方程的系數(shù)矩陣A和常數(shù)項(xiàng)L，組成法方程式。 （8）解法方程，求得外方位元素的改正數(shù)dXs，dYs，dZs，dq，dw，dk。 （9）用前次迭代取得的近似值，加本次迭代的改正數(shù)，計(jì)算外方位元素的新值。（10） 將求得的外方位元素改正數(shù)與規(guī)定的限差比較，若小于限差則迭代結(jié)束。否則用新的近似值重復(fù)（4）（9），直到滿足要求為止。四程序流程圖。 五程序的主要源代碼如下所示： #include "stdafx.h"#include "KongJianHouFangJiao

5、Hui.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()/導(dǎo)入測量數(shù)據(jù)按鈕代碼CFileDialog

6、dlgOpenFile(TRUE, _T("txt"), NULL,OFN_FILEMUSTEXIST, _T("(文本文件)|*.txt|(所有文件）|*.*)|");if (dlgOpenFile.DoModal() = IDCANCEL)return;/如果選擇取消按鈕，則退出CString strPathName = dlgOpenFile.GetPathName();/獲取選擇的文件的完整路徑CString strFileTitle = dlgOpenFile.GetFileTitle();/獲取文件名CString strExtName =

7、dlgOpenFile.GetFileExt();/獲取文件擴(kuò)展名CStdioFile sf; /創(chuàng)建文件對(duì)象/以讀的形式打開文件，如果打開失敗if (!sf.Open(strPathName, CFile:modeRead)MessageBox(_T("讀取文件出錯(cuò)！");return;CString strLine;/存放每一行文本CStringArray strSplit;CeHuiLei SplitString;int hang = 0;while (sf.ReadString(strLine)if (hang <= 1)hang+;else int n =

8、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()

9、;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;/表示導(dǎo)入數(shù)據(jù)按鈕已被點(diǎn)擊并且已導(dǎo)入數(shù)據(jù)CMatrix CKongJianHouFangJiaoHuiDlg:GetX(CMatrix A, CMatrix L)/計(jì)算外方位元素改正數(shù)矩陣CM

10、atrix X,_A,_AA,N_AA;_A = A;/A的轉(zhuǎn)置_AA = _A*A;N_AA = _AA.Inv();/_AA的逆矩陣X = N_AA*_A*L;return X;CMatrix CKongJianHouFangJiaoHuiDlg:GetA(CMatrix xyXYZ, double f, CMatrix XX)/計(jì)算系數(shù)矩陣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,

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

12、 = -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) +

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

14、, 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

15、)/計(jì)算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

16、);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

17、 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*(

18、YA - YS) + c3*(ZA - ZS);L(2*i, 0) =x - x0;L(2*i+1, 0) =y - y0;return L;void CKongJianHouFangJiaoHuiDlg:OnBnClickedButton1()/計(jì)算按鈕UpdateData(true);if (ButtonClickeOk = 1)MessageBox(_T("請(qǐng)先點(diǎn)擊“導(dǎo)入已知坐標(biāo)”按鈕導(dǎo)入數(shù)據(jù)！"), _T("提示");return;else if (strf = _T("") | strm = _T("")

19、 | strx0 = _T("") | stry0 = _T("")MessageBox(_T("請(qǐng)輸入數(shù)據(jù)！"), _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

20、);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);/迭代計(jì)算：do CMatrix A = GetA(xyXYZ, _wtof(strf), XX);CMatrix L = GetL(xyXYZ, _wtof(strf), XX);dX = GetX(A, L);for (int n = 0; n < 6; n+)XX(0, n) = XX(0, n) + dX(n, 0); while (dX(3, 0) >= 2.908882087e-5 | dX(4, 0) >= 2.908882087e-5 | dX(5, 0) >= 2.908882087e-5);/各角元素迭代計(jì)算至其改正