摄影测量学空间后交前交实验报告

1、中南大学本科生课程设计（实践）任务书、设计报告（摄影测量与遥感概论）题 目： 空间后方交会-前交 院 系：地球科学与信息物理学院班 级： 测绘1201班 指 导 老 师：邹峥嵘 学 号： 0103120506 姓 名： 杨腊梅 二零一四年十一月一、实验目的通过对数字影像空间后交前交的程序设计实验，要求我们进一步理解和掌握影像外方位元素的有关理论、原理和方法。利用计算机程序设计语言编写摄影测量空间交会软件进行快速确定影响的外方位元素及其精度，然后通过求得的外方位元素求解位置点的地面摄影测量坐标，达到通过摄影测量量测地面地理坐标的目的。二、实验要求1.用C、VB、C+或MATLAB语言编写空间后方

2、交会-空间前方交会程序2.提交实习报告：程序框图、程序源代码、计算结果、体会3.计算结果：像点坐标、地面坐标、单位权中误差、外方位元素及其精度4.完成时间：11月11日前完成三、实验数据点号左片右片地面摄影测量坐标系xyxyXYZGCP116.01279.963-73.9378.7065083.2055852.099527.925GCP288.5681.134-5.25278.1845780.025906.365571.549GCP313.362-79.37-79.122-78.8795210.8794258.446461.81GCP482.24-80.027-9.887-80.0895909

3、.2644314.283455.484151.75880.555-39.95378.463214.618-0.231-76.0060.036349.88-0.782-42.201-1.022486.14-1.346-7.706-2.112548.035-79.962-44.438-79.736f=150.000mm，x0=0，y0=0四、实验思路Ø 利用后方交会得出两张相片各自的外方位元素1）获取已知数据影响比例尺m,，内方位元素x0 、y0 、f ，控制点的地面摄影测量坐标Xtp, Ytp, Ztp2）量测控制点左片和右片的像点坐标 x，y3）确定未知数初值 Xs0, Ys0, Z

4、s0, ，（线元素可用控制点均值代替，角元素可用0初始化），即：，=0，=0，=04）计算旋转矩阵R5）利用共线方程逐点计算像点坐标的近似值6）组成误差方程式并法化7）解求外方位元素改正数8）检查迭代是否收敛（改正值是否小于某一特定常数）Ø 利用解出的外方位元素进行前方交会1）获取已知数据x0 , y0 , f, XS1, YS1, ZS1,1,1,1 , XS2, YS2, ZS2，2，2，2 2）量测像点坐标 x1,y1 ,x2,y2 3）由外方位线元素计算基线分量BX, BY, BZ 4）由外方位角元素计算像空间辅助坐标 X1, Y1, Z1 , X2, Y2, Z2 5）计算

5、点投影系数 N1 , N2 6）计算地面坐标 XA, YA, ZA 五、实验过程Ø 程序流程图此过程完成空间后方交会求解像片的外方位元素，其中改正数小于限差（0.00003，相当于0.1的角度值）为止。在这个过程中采用迭代的方法，是外方位元素逐渐收敛于理论值，每次迭代所得的改正数都应加到上一次的初始值之中。Ø 程序中的类与其主要成员、主要函数设计ControlPointData类：控制点、未知点的点名、地面摄影测量坐标、左右像片坐标Photo类：左右片的内方位元素、外方位元素及其精度Fileheader类：相片对象、控制点个数、未知点个数、控制点数组、未知点数组Adjust

6、ment类：头文件对象、控制点数组、未知点数组、控制点个数、未知点个数、 bool LoadData(const CString& strFileName,CString& strData);/读取已知函数CString* SplitString(CString str,char split, int& iSubStrs);/分割函数void adjustment1(CString &OutData);/后方交会函数void adjustment2(CString &OutData);/前方交会函数void SaveData(CString &O

7、utData);/保存结果Ø 读取数据文件格式（.txt）elements of interior orientation,0,0,150.000GivenControlPoint,4GCP1,16.012,79.963,-73.93,78.706,5083.205,5852.099,527.925 GCP2,88.56,81.134,-5.252,78.184,5780.02,5906.365,571.549 GCP3,13.362,-79.37,-79.122,-78.879,5210.879,4258.446,461.81 GCP4,82.24,-80.027,-9.887,-

8、80.089,5909.264,4314.283,455.484 UnKnownPoint,51,51.758,80.555,-39.953,78.463 2,14.618,-0.231,-76.006,0.036 3,49.88,-0.782,-42.201,-1.022 4,86.14,-1.346,-7.706,-2.112 5,48.035,-79.962,-44.438,-79.736 End DataØ 程序代码ControlPointData.hclass ControlPointData : public CDialogpublic:CString strID;/控制

9、点名double x1,y1,x2,y2;/左片、右片像点坐标 double Xtp,Ytp,Ztp;/地面摄影测量坐标double Xs,Ys,Zs;/物方空间坐标double ;/单位权中误差 bool state;class Photospublic:int x0,y0,f,m;/内方位元素 double ,Xs1,Xs2,Ys1,Ys2,Zs1,Zs2;/外方位元素double ;/单位权中误差 double Xs1,Xs2,Ys1,Ys2,Zs1,Zs2;double ,;/协因数Photos(void)x0=0;y0=0;f=0;class fileheader public: P

10、hotos photo;/内方位元素int StationCount;/控制点个数int OriAngleCount;/测站个数 int UnKnownPointCount;/待定点个数ControlPointData *StationCoor;/控制点坐标信息ControlPointData *StationCoorJ;/控制点坐标信息ControlPointData *UnKnownPoint;/待定点坐标信息; Adjustment.h#pragma once#include "ontrolPointData.h"#include <locale.h> #

11、include"math.h" #include "Matrix.h"class Adjustmentpublic:fileheader header;/头文件 ControlPointData *KnownPoint;/已知点int KnownPointCount;/已知点个数int UnKnownPointCount;/待定点个数public:Adjustment(void);Adjustment(void);public:bool LoadData(const CString& strFileName,CString& strDat

12、a);/读取已知函数CString* SplitString(CString str,char split, int& iSubStrs);/分割函数public:void adjustment1(CString &OutData);/后方交会函数void adjustment2(CString &OutData);/前方交会函数void SaveData(CString &OutData);/保存结果;Adjustment.cpp#include "StdAfx.h"#include "Adjustment.h"#inc

13、lude "ontrolPointData.h"#include "Space resection and forward intersection.h"#include "Space resection and forward intersectionDlg.h"#include "Adjustment.h"#include <locale.h> #include"math.h" #include "Matrix.h"#define pai 3.1415926C

14、String* Adjustment:SplitString(CString str,char split, int& iSubStrs)/分割函数 int iPos = 0; /分割符位置 int iNums = 0; /分割符的总数 CString strTemp = str; CString strRight; /先计算子字符串的数量 while (iPos != -1) iPos = strTemp.Find(split); if (iPos = -1) break; strRight = strTemp.Mid(iPos + 1, str.GetLength(); strTe

15、mp = strRight; iNums+; if (iNums = 0) /没有找到分割符 /子字符串数就是字符串本身 iSubStrs = 1; return NULL; /子字符串数组 iSubStrs = iNums + 1; /子串的数量= 分割符数量+ 1 CString* pStrSplit; pStrSplit = new CStringiSubStrs; strTemp = str; CString strLeft; for (int i = 0; i < iNums; i+) iPos = strTemp.Find(split); /左子串 strLeft = str

16、Temp.Left(iPos); /右子串 strRight = strTemp.Mid(iPos + 1, strTemp.GetLength(); strTemp = strRight; pStrSpliti = strLeft; pStrSplitiNums = strTemp; return pStrSplit;bool Adjustment:LoadData(const CString& strFileName,CString& strData)CStdioFile sf; if(!sf.Open(strFileName, CFile:modeRead) return

17、 false; CString strLine; int n1; CString *strTmp=NULL; while(strLine!=_T("End Data")/读取文件头 sf.ReadString (strLine); strData+=strLine; strData+=_T("rn");if(strLine!=_T("End Data") strTmp = SplitString(strLine, ',',n1); else break; if(strTmp0="elements of int

18、erior orientation") header.photo.x0=_tstof(strTmp1);header.photo.y0=_tstof(strTmp2);header.photo.f=_tstof(strTmp3);else if(strTmp0="GivenControlPoint") header.StationCount =_ttoi(strTmp1);header.StationCoor =new ControlPointDataheader.StationCount; if(strTmp!=NULL)/释放内存 delete strTmp;

19、 strTmp=NULL; for(int i=0;i<header.StationCount;i+) sf.ReadString (strLine); strData+=strLine; strData+=_T("rn"); strTmp = SplitString(strLine, ',',n1); header.StationCoori.strID=strTmp0;header.StationCoori.x1 =_tstof(strTmp1);header.StationCoori.y1=_tstof(strTmp2);header.Statio

20、nCoori.x2 =_tstof(strTmp3);header.StationCoori.y2 =_tstof(strTmp4);header.StationCoori.Xtp=_tstof(strTmp5);header.StationCoori.Ytp =_tstof(strTmp6);header.StationCoori.Ztp=_tstof(strTmp7); if(strTmp!=NULL)/释放内存 delete strTmp; strTmp=NULL; else if(strTmp0="UnKnownPoint") header.UnKnownPoint

21、Count =_ttoi(strTmp1);header.UnKnownPoint=new ControlPointDataheader.UnKnownPointCount; if(strTmp!=NULL)/释放内存 delete strTmp; strTmp=NULL; for(int i=0;i<header.UnKnownPointCount;i+) sf.ReadString (strLine); strData+=strLine; strData+=_T("rn"); strTmp = SplitString(strLine, ',',n1

22、); header.UnKnownPointi.strID=strTmp0;header.UnKnownPointi.x1 =_tstof(strTmp1);header.UnKnownPointi.y1=_tstof(strTmp2);header.UnKnownPointi.x2 =_tstof(strTmp3);header.UnKnownPointi.y2 =_tstof(strTmp4);header.UnKnownPointi.Xtp=0;header.UnKnownPointi.Ytp =0;header.UnKnownPointi.Ztp=0; if(strTmp!=NULL)

23、/释放内存 delete strTmp; strTmp=NULL; if(strTmp!=NULL)/释放内存delete strTmp;strTmp=NULL; sf.Close();/关闭文件return true;const double pi=3.1415926; void Adjustment:adjustment1(CString &OutData)/单片后交 CMatrix B,L,X;B.SetSize(2*header.StationCount,6);/系数阵 L.SetSize(2*header.StationCount,1);/常数项阵X.SetSize(6,1)

24、;/未知数改正数之和的阵CMatrix R(3,3),Xb(3,1),bb(3,1);/R为旋转矩阵，Xb为X-Xs,Y-Ys,Z-Zs的矩阵，bb为X-，Y-,Z-的矩阵CMatrix Dx(6,1),Nbb,Nvv;/Dx为改正数阵header.photo.m=10000;double Xs0,Ys0,Zs0,;/未知数初值/确定未知数初值double x=0,y=0,z=0;for(int i=0;i<header.StationCount;i+)x+=header.StationCoori.Xtp;y+=header.StationCoori.Ytp;z+=header.Stat

25、ionCoori.Ztp;Xs0=x/header.StationCount;Ys0=y/header.StationCount;Zs0=header.photo.m*header.photo.f/1000;=0;=0;=0;double n=0;/计算左片外方位元素*/*/don+;/计算旋转矩阵Rdouble a1,a2,a3,b1,b2,b3,c1,c2,c3;a1=cos()*cos()-sin()*sin()*sin(); R(0,0)=a1;a2=-cos()*sin()-sin()*sin()*cos(); R(0,1)=a2;a3=-sin()*cos(); R(0,2)=a3

26、;b1=cos()*sin(); R(1,0)=b1;b2=cos()*cos(); R(1,1)=b2;b3=-sin(); R(1,2)=b3;c1=sin()*cos()+cos()*sin()*sin(); R(2,0)=c1;c2=-sin()*sin()+cos()*sin()*cos(); R(2,1)=c2;c3=cos()*cos(); R(2,2)=c3;/计算像点坐标的近似值,构建B，L矩阵for(int i=0; i<header.StationCount;i+)bb(0,0)=a1*(header.StationCoori.Xtp-Xs0)+b1*(header

27、.StationCoori.Ytp-Ys0)+c1*(header.StationCoori.Ztp-Zs0);bb(1,0)=a2*(header.StationCoori.Xtp-Xs0)+b2*(header.StationCoori.Ytp-Ys0) /单位为m+c2*(header.StationCoori.Ztp-Zs0);bb(2,0)=a3*(header.StationCoori.Xtp-Xs0)+b3*(header.StationCoori.Ytp-Ys0)+c3*(header.StationCoori.Ztp-Zs0); double t1=header.Statio

28、nCoori.x1;/（X-X0）(Y-Y0)的值，单位为mmdouble t2=header.StationCoori.y1;L(2*i,0)=(header.StationCoori.x1+header.photo.f*bb(0,0)/bb(2,0)/1000;L(2*i+1,0)=(header.StationCoori.y1+header.photo.f*bb(1,0)/bb(2,0)/1000;B(2*i,0)=(a1*header.photo.f+a3*t1)/bb(2,0)/1000;B(2*i,1)=(b1*header.photo.f+b3*t1)/bb(2,0)/1000;

29、B(2*i,2)=(c1*header.photo.f+c3*t1)/bb(2,0)/1000;B(2*i,3)=(t2*sin()-(t1/header.photo.f*(t1*cos()-t2*sin()+header.photo.f*cos()*cos()/1000;B(2*i,4)=(-header.photo.f*sin()-t1/header.photo.f*(t1*sin()+t2*cos()/1000;B(2*i,5)=t2/1000;B(2*i+1,0)=(a2*header.photo.f+a3*t2)/bb(2,0)/1000;B(2*i+1,1)=(b2*header.

30、photo.f+b3*t2)/bb(2,0)/1000;B(2*i+1,2)=(c2*header.photo.f+c3*t2)/bb(2,0)/1000;B(2*i+1,3)=(-t1*sin()-(t2/header.photo.f*(t1*cos()-t2*sin()-header.photo.f*sin()*cos()/1000;B(2*i+1,4)=(-header.photo.f*cos()-t2/header.photo.f*(t1*sin()+t2*cos()/1000;B(2*i+1,5)=-t1/1000;/平差计算Nbb=(B)*B;Nvv=(B)*L;Dx=(Nbb.I

31、nv()*Nvv;Xs0+=Dx(0,0);Ys0+=Dx(1,0);Zs0+=Dx(2,0);+=Dx(3,0);+=Dx(4,0);+=Dx(5,0);while(fabs(Dx(0,0)>0.001)|(fabs(Dx(1,0)>0.001)|(fabs(Dx(2,0)>0.001)|(fabs(Dx(3,0)*206265)>1)|(fabs(Dx(4,0)*206265)>1)|(fabs(Dx(5,0)*206265)>1);header.photo.Xs1=Xs0;header.photo.Ys1=Ys0;header.photo.Zs1=Zs

32、0;header.photo.=;header.photo.=;header.photo.=;CMatrix V=B*Dx-L;CMatrix Vt=(V)*V;/改正数转置 double m0=sqrt(Vt(0,0)/(2*header.StationCount-6); CMatrix Q=Nbb.Inv();/协因数阵 /求解点位误差 header.photo.Xs1=sqrt(Q(0,0)*m0; header.photo.Ys1=sqrt(Q(1,1)*m0; header.photo.Zs1=sqrt(Q(2,2)*m0;header.photo.=sqrt(Q(3,3)*m0;

33、header.photo.=sqrt(Q(4,4)*m0; header.photo.=sqrt(Q(5,5)*m0; CString data;data.Format(_T("%srn"),_T("后方交会成果");OutData+=data;OutData+="rn"data.Format(_T("%s %s %.f"),_T("左片"),_T("迭代次数"),n);OutData+=data; OutData+=("rnrn");data.Forma

34、t(_T(" %srn %s %s %s %s %s %srn"),_T("左片外方位元素(线元素为m,角元素为度)"),_T("Xs"),_T("Ys"),_T("Zs"),_T(""),_T(""),_T(""); OutData+=data; data.Format(_T("%.4f %.4f %.4f %.4f %.4f %.4frnrn"),header.photo.Xs1,header.photo.Ys1

35、,header.photo.Zs1,header.photo.*180/pai,header.photo.*180/pai,header.photo.*180/pai);OutData+=data;data.Format(_T("%srn%s %.3frnrn %s %s %s %s %s %srn"),_T("左片精度(线元素精度为m,角元素精度为）"),_T("单位权中误差m0"),m0,_T("Xs"),_T("Ys"),_T("Zs"),_T(""

36、),_T(""),_T(""); OutData+=data; data.Format(_T(" %.4f %.4f %.4f %.4f %.4f %.4frn"), header.photo.Xs1,header.photo.Ys1,header.photo.Zs1,header.photo.*180/pai*3600,header.photo.*206265,header.photo.*206265);OutData+=data; OutData+=("rnrn");/计算右片外方位元素*/*/Xs0=x/he

37、ader.StationCount;Ys0=y/header.StationCount;Zs0=header.photo.m*header.photo.f/1000;=0;=0;=0;n=0;don+;/计算旋转矩阵Rdouble a1,a2,a3,b1,b2,b3,c1,c2,c3;a1=cos()*cos()-sin()*sin()*sin(); R(0,0)=a1;a2=-cos()*sin()-sin()*sin()*cos(); R(0,1)=a2;a3=-sin()*cos(); R(0,2)=a3;b1=cos()*sin(); R(1,0)=b1;b2=cos()*cos();

38、 R(1,1)=b2;b3=-sin(); R(1,2)=b3;c1=sin()*cos()+cos()*sin()*sin(); R(2,0)=c1;c2=-sin()*sin()+cos()*sin()*cos(); R(2,1)=c2;c3=cos()*cos(); R(2,2)=c3;/计算像点坐标的近似值,构建B，L矩阵for(int i=0; i<header.StationCount;i+)bb(0,0)=a1*(header.StationCoori.Xtp-Xs0)+b1*(header.StationCoori.Ytp-Ys0)+c1*(header.StationC

39、oori.Ztp-Zs0);bb(1,0)=a2*(header.StationCoori.Xtp-Xs0)+b2*(header.StationCoori.Ytp-Ys0) /单位为m+c2*(header.StationCoori.Ztp-Zs0);bb(2,0)=a3*(header.StationCoori.Xtp-Xs0)+b3*(header.StationCoori.Ytp-Ys0)+c3*(header.StationCoori.Ztp-Zs0); double t1=header.StationCoori.x2;/（X-X0）(Y-Y0)的值，单位为mmdouble t2=h

40、eader.StationCoori.y2;L(2*i,0)=(header.StationCoori.x2+header.photo.f*bb(0,0)/bb(2,0)/1000;L(2*i+1,0)=(header.StationCoori.y2+header.photo.f*bb(1,0)/bb(2,0)/1000;B(2*i,0)=(a1*header.photo.f+a3*t1)/bb(2,0)/1000;B(2*i,1)=(b1*header.photo.f+b3*t1)/bb(2,0)/1000;B(2*i,2)=(c1*header.photo.f+c3*t1)/bb(2,0)

41、/1000;B(2*i,3)=(t2*sin()-(t1/header.photo.f*(t1*cos()-t2*sin()+header.photo.f*cos()*cos()/1000;B(2*i,4)=(-header.photo.f*sin()-t1/header.photo.f*(t1*sin()+t2*cos()/1000;B(2*i,5)=t2/1000;B(2*i+1,0)=(a2*header.photo.f+a3*t2)/bb(2,0)/1000;B(2*i+1,1)=(b2*header.photo.f+b3*t2)/bb(2,0)/1000;B(2*i+1,2)=(c2

42、*header.photo.f+c3*t2)/bb(2,0)/1000;B(2*i+1,3)=(-t1*sin()-(t2/header.photo.f*(t1*cos()-t2*sin()-header.photo.f*sin()*cos()/1000;B(2*i+1,4)=(-header.photo.f*cos()-t2/header.photo.f*(t1*sin()+t2*cos()/1000;B(2*i+1,5)=-t1/1000;/平差计算Nbb=(B)*B;Nvv=(B)*L;Dx=(Nbb.Inv()*Nvv;Xs0+=Dx(0,0);Ys0+=Dx(1,0);Zs0+=Dx

43、(2,0);+=Dx(3,0);+=Dx(4,0);+=Dx(5,0);while(fabs(Dx(0,0)>0.001)|(fabs(Dx(1,0)>0.001)|(fabs(Dx(2,0)>0.001)|(fabs(Dx(3,0)*206265)>1)|(fabs(Dx(4,0)*206265)>1)|(fabs(Dx(5,0)*206265)>1);header.photo.Xs2=Xs0;header.photo.Ys2=Ys0;header.photo.Zs2=Zs0;header.photo.=;header.photo.=;header.pho

44、to.=;V=B*Dx-L; Vt=(V)*V;/改正数转置 m0=sqrt(Vt(0,0)/(2*header.StationCount-6); Q=Nbb.Inv();/协因数阵 /求解点位误差 header.photo.Xs2=sqrt(Q(0,0)*m0; header.photo.Ys2=sqrt(Q(1,1)*m0; header.photo.Zs2=sqrt(Q(2,2)*m0;header.photo.=sqrt(Q(3,3)*m0; header.photo.=sqrt(Q(4,4)*m0; header.photo.=sqrt(Q(5,5)*m0; data.Format(

45、_T("%s %s %.f"),_T("右片"),_T("迭代次数"),n);OutData+=data; OutData+=("rnrn");data.Format(_T(" %srn %s %s %s %s %s %srn"),_T("右片外方位元素(线元素为m,角元素为度)"),_T("Xs"),_T("Ys"),_T("Zs"),_T(""),_T(""),_T(&qu

46、ot;"); OutData+=data; data.Format(_T("%.4f %.4f %.4f %.4f %.4f %.4frnrn"),header.photo.Xs2,header.photo.Ys2,header.photo.Zs2,header.photo.*180/pai,header.photo.*180/pai,header.photo.*180/pai);OutData+=data;data.Format(_T("%srn%s %.4f"),_T("右片精度(线元素精度为m,角元素精度为)"),_T

47、("单位权中误差m0（单位为um)"),m0);OutData+=data; OutData+=("rnrn");data.Format(_T(" %s %s %s %s %s %srn"),_T("Xs"),_T("Ys"),_T("Zs"),_T(""),_T(""),_T(""); OutData+=data; data.Format(_T(" %.4f %.4f %.4f %.4f %.4f %.4frnrnrn"), header.photo.Xs2,header.photo.Ys2,header.photo.Zs2,header.photo.*180/pai*3600,header.photo.*206265,header.photo.*206265);OutData+=