测绘程序设计基础课程设计报告

1、中南大学测绘程序设计基础课程设计报告专业班级测绘1001班学生姓名李佳霖学生学号0405100111学 院地球科学与信息物理学院指导教师范冲老师2012年7月30日实验目的使学生能够在学习完测绘程序程序设计基础（C+）、测量学、测量平差等课程的基础上，进一步巩固这些知识，并培养学生的综合应用能力。实验内容及要求 内容：设计一个导线网平差的程序。.要求：课程设计的第一、二天根据课堂上所讲进行系统设计，对主要的概算与平差计算要求列出计算的详细步骤，并写出流程图，最后形成系统设计说明书，并以天为单位制定详细的程序设计计划。进行编程，要求编程规范编写代码。 = 3 * GB3 考核：要求讲述编程思路与

2、具体代码的功能含义，并演示程序。系统设计1、需求分析控制网平差程序对野外控制网观测数据进行平差数据处理，其目的就是根据最小二乘原理，消除网中的各种几何矛盾，求出全网各待定元素（未知点的平面坐标或三维坐标）。2、实现功能 1).数据输入：使用外部文件的导入，读取.txt格式的数据文件实现数据的读取。 2).概算：主要实现近似坐标的概算。 3).平差计算：对观测数据进行精密平差计算，得到平差后的点位坐标，方向观测值，边长观测值等，精度评定。 4).成果输出：控制网图形输出，平差结果报表及其输出，绘制误差椭圆等。3、系统总体框架边角网平差程序近似坐标概算成果输出平差计算数据输入四、实现思路1、数据输

3、入模块、数据格式2已知数据， 2为已知点的个数 点名,X,Y 未知点数据，3为未知点的个数 点名33,4,58观测数据18为方向观测值的个数；8代表距离观测值个数。18.、数据输入流图数据导入描述的.txt格式数据文件数据编辑数据输入外业数据采集、具体实现 先是定义一个控制点类CControlPoint，用来定义控制点m_pKnownPoint/m_pUnknownPoint、方向角观测值CBearingObs、距离观测值CDistObs的格式，这个类编译好之后，就可以用于数据的读取了。CControlPoint* m_pUnknownPoint;/未知点数组int m_iUnknownPoi

4、ntCount; /未知点个数CControlPoint* m_pKnownPoint; /已知点数组int m_iKnownPointCount; /已知点个数CBearingObs* m_pBearingObs; /方向观测值数组int m_iBearingObsCount; /方向观测值个数CDistObs* m_pDistObs; /距离观测值数组int m_iDistObsCount; /距离观测值个数CBearingObs* m_pZeroBearing; /零方向观测值数据int m_iZeroBearingCount; /零方向观测值个数bool LoadData(const

5、CString& strFileName);CString* SplitString(CString str , char split, int iSubStrs);/字符串分割函数按规定格式编辑好数据文件后，将这个txt文档导入，通过MFC文件读取方法逐行读入数据，然后通过调用分割函数SplitString，逐行进行分割，将数据保存到定义的动态数组中，这样就实现了数据的读取与保存。2、概算模块（主要是近似坐标的推算）近似坐标推算公式 EQ oac(,1)：按方向观测值计算三角网中待定点的近似坐标 EQ oac(,2)按边长观测值计算三边网中待定点的近似坐标先按公式(2)计算三角形的内角，再用

6、公式（1）计算待定点的坐标 EQ oac(,3)按方向和边长观测值计算导线网中待定点的近似坐标 综合比较上述三种推算近似坐标的方法，第一和第二种方法，在近似坐标推算时必须有三角形才能计算，适用于三角网。而本次试验采用的是边角网的观测数据，对于第一种只按角推算和第二种只按边推算，显然也是不合适的。而第三种极坐标推算的方法，对所有网形均适用，从通用性和数据合理利用这两面考虑，故采用第三种极坐标推算的方法推算近似坐标。 在具体计算中采用逐点解算法，选择已知点或已计算出坐标的点作为起算点，根据具体的图形，来逐个推算出各点的坐标。这种方法比较容易实现，而对于整体解算的方法，对所有未知点进行整体解算，实现

7、起来太困难。推算流程图如下：3、平差计算模块 EQ oac(,1)附有条件的间接平差原理：设未知数个数为 t,观测值个数为n,条件方程个数为m误差方程：V=BX+L条件方程：AX+W=0根据最小二乘原理：在条件AX +W=0下，即求条件极值法，利用拉格朗日不定乘数法：由于此次实验中间接平差更适合编程的特点，故我们选用间接平差。控制网平差程序对野外控制网观测数据进行平差数据处理，其目的就是根据最小二乘原理，消除网中的各种几何矛盾，求出全网各待定元素（未知点的平面坐标或三维坐标）。 EQ oac(,2)平差计算步骤：读取观测数据与已知点坐标数据；（2）计算未知点近似坐标（用“导线点坐标增量法近似坐

8、标”算法）；（3）计算角度观测值误差方程的系数和常数项；（4）计算边长观测值误差方程的系数和常数项；（5）列条件式；（6）组成误差方程，得到系数矩阵和常数项矩阵；（7）定权，得到权矩阵；（8）组成方程；（9）解算法方程，求得dX=dx1 dy1 dx2 dy2 （10）求得平差后的坐标值X=X+dX。因为只有dx 很小时误差方程式、条件方程式才是严密的，因此当dX 的值较大时应把X=X_+dX 作为新的近似值重复1-3 步，直到dx 中绝对值最大的一个小于给定的限值。（11）精度评定，计算误差椭圆参数等 EQ oac(,3)误差方程式及条件方程式的列立以未知点坐标为平差对象，一个观测值对应一个

9、误差方程式，误差方程式包括两类：方向观测误差方程式和边长观测误差方程式。方向观测误差方程式 （2）边长观测误差方程式 EQ oac(,4)定权 根据先验方向观测中误差和边长观测中误差来定权（1） 方向观测值的权 在等精度方向观测的控制网中，可把方向观测值的权定为1 （2）边长观测值的权而Ms=A+B*S/1000;其中S单位为Km,这样的话Ms单位为mm EQ oac(,5)平差后的坐标值计算求得平差后的坐标值X=X+dX。因为只有dx很小时误差方程式、条件方程式才是严密的，因此当dX的值较大时应把X=X_+dX作为新的近似值进行迭代，直到dx中绝对值最大的一个小于给定的0.1mm（注意：这里

10、只是指坐标改正值，不包括定向改正值，所以在求dx中绝对值最大值时，应除去定向改正值，否则会无限迭代计算，因为定向改正值很大，而且迭代计算时不会进行改正）。 EQ oac(,6)精度评定中误差计算：，n为观测值个数，t为未知数个数待定点点位误差:误差椭圆元素计算:长半径方位角长半径、短半径数据读取，并保存到相应数组中平差计算流程图调用坐标概算函数，推算待求点近似坐标调用调用计算边长观测值误差方程的系数和常数项及权计算角度观测值误差方程的系数和常数项及权组成误差方程，得到系数矩阵和常数项矩阵进行间接平差否If(MaxAbsDetX(detx)0.1)是精度评定平差成果输出结束成果输出模块成果输出

11、控制网图形误差椭圆的绘制控制点成果精度统计具体实现：平差的成果按一定的格式格式化输出到界面上，并以txt文档的格式保存。调用绘制控制网图函数和绘制误差椭圆函数实现图形的输出显示，新建一Dlg用于设置误差椭圆的比例系数。五、界面设计图中主要有两个按钮，Open the file and Compute按钮和Draw the Ellipse按钮分别用作打开文件并计算输出和绘制误差椭圆。六、实验总结通过本次实验我在学习完测绘程序程序设计基础（C+）、测量学、测量平差等课程的基础上，通过实践巩固这些课程上的理论知识，并培养了动手能力，综合应用的能力。在本次实验当中，主要设计思路参考了之前上课时所编写的

12、水准网平差和单导线平差的思路，其中数据读入模块方法基本一致，需要编写坐标概算和平差计算以及控制网图形及误差椭圆的绘制工作。在实验中遇到了一下的几个问题：1、坐标概算考虑不周，在坐标概算时，只考虑到了已知数据在前面给出的情况，而如果已知数据出现在后面则会出现bug，在范老师的指导下，发现并改正了这个小问题；2、平差计算中出现问题，在平差计算时发现改正以后与正确值相差更多了，我自己检查发现在组成法方程时，由于自己疏忽，写错了代码，改过之后就能够输出了正确结果；3、绘制控制网图形和误差椭圆时，以前的方法是建立单文档程序，在窗口中画图，但是我这次这样做发现在这个里边以前程序中的部分功能实现不了，我想应

13、该是因为那些函数是基于对话框的，所以就将画图的功能放入对话框中实现了。 衷心感谢范老师张老师和学长在实验当中对我的指导。附：主要代码#pragma once#include ControlPoint.h#include Angle.h#include Matrix.hconst double Ro = 206264.806247; /一弧度对应的秒数class CTraverseNetAdjpublic:CTraverseNetAdj(void);CTraverseNetAdj(void);private: CControlPoint* m_pUnknownPoint;/未知点数组int m_

14、iUnknownPointCount; /未知点个数 CControlPoint* m_pKnownPoint; /已知点数组int m_iKnownPointCount; /已知点个数CBearingObs* m_pBearingObs; /方向观测值数组int m_iBearingObsCount; /方向观测值个数 CDistObs* m_pDistObs; /距离观测值数组int m_iDistObsCount; /距离观测值个数CBearingObs* m_pZeroBearing; /零方向观测值数据int m_iZeroBearingCount; /零方向观测值个数CString

15、* SplitString(CString str , char split, int &iSubStrs);/分割字符串void ApproximateXY();/坐标概算，求近似坐标CControlPoint* SearchKnownPointUsingNAME(CString strNAME);/根据点名查找已知点，返回该点的指针CControlPoint* SearchUnknownPointUsingNAME(CString strNAME);/根据点名查找未知知点，返回该点的指针CControlPoint* SearchPointUsingNAME(CString strNAME)

16、;/根据点号查找控制点（包括已知点和未知点），返回该点的指针public:void SetUnknownPointSize(int size);/设置未知点数组大小void SetKnownPointSize(int size);/设置已知点数据大小void SetBearingObsSize(int size); /设置方向观测值数组大小void SetDistObsSize(int size); /设置距离观测值数组大小void SetZeroBearingSize(int size);/设置零方向观测值数组大小CAngle ZeroBearing(const CBearingObs m_

17、BearingObs); /计算定向零度方向角函数输入一个方向观察值返回测站点的零度定向角public:CAngle Azi(const CControlPoint& P1, const CControlPoint& P2);/已知两点求方位角，输入参数为CControlPoint类型double HDist(const CControlPoint& P1, const CControlPoint& P2);/求两点之间的距离bool LoadObsData(const CString& strFileName);/从文件中导入数据 void RigorousAdjust(const CStr

18、ing& strFileName);/严密平差（间接平差）计算主函数 private:double MaxAbsDetX(const CMatrix& detX);/求未知改正数绝对值最大的void OutMatrixToFile(const CMatrix& mat,CStdioFile& SF);/把矩阵输出到文件中CMatrix Weight(double SigmaAng, double SigmaDist);/简单定权 /计算方向观测方程的系数和常数项void ComBearingObsEff(const CBearingObs& BearingObs, double& Li,dou

19、ble& Aj,double& Bj, double& Ak,double& Bk); /计算距离观测方程的系数和常数项void ComDistObsEff(const CDistObs& DistObs, double& Li, double& Ci, double& Di, double& Ck, double& Dk); /输入一个方向观测对象，查找这一方向观测对象所对应的零方向对象在零方向数组中的位置int SearchZeroBearing(const CBearingObs& BearingObs); /组成误差方程,B 为系数矩阵，L为常数项向量 void FormErrorEq

20、uations(CMatrix& B, CMatrix& L);/画图public:double dScale; /比例/坐标原点double dOrgX; double dOrgY;int iLineWidth;/线宽COLORREF crColor ;public:void SetScale(double scale);/设置绘图比例 void SetOrg(double x, double y);/设置坐标原点 void SetLineWidth(int iWidth);/设置线宽 void SetColor(COLORREF color);/设置颜色void draw(CDC* pDC

21、, CRect& rect);/绘制图形;#include StdAfx.h#include TraverseNetAdj.h#include CommonSurveyFunctions.h#include math.h#include CTraverseNetAdj:CTraverseNetAdj(void)/释放动态数组内存if(m_pUnknownPoint!=NULL)delete m_pUnknownPoint;m_pUnknownPoint=NULL;if(m_pKnownPoint!=NULL)delete m_pKnownPoint;m_pKnownPoint=NULL; if

22、(m_pBearingObs!=NULL)delete m_pBearingObs;m_pBearingObs=NULL; if(m_pDistObs!=NULL)delete m_pDistObs;m_pDistObs=NULL;if (m_pZeroBearing!=NULL)delete m_pZeroBearing;m_pZeroBearing=NULL;/导线网的已知点个数不确定，所以不能像单导线那样直接确定已知点个数m_pKnownPoint=NULL; /已知点数组 m_iKnownPointCount=0; /已知点个数m_pUnknownPoint=NULL; /未知点数组m

23、_iUnknownPointCount=0; /未知点个数 m_pBearingObs=NULL; /角度观测值数组m_iBearingObsCount=0; /角度观测值个数 m_pDistObs=NULL; /距离观测值数组m_iDistObsCount=0; /距离观测值个数m_pZeroBearing=NULL; /零方向观测值数组m_iZeroBearingCount=0; /零方向观测值个数CTraverseNetAdj:CTraverseNetAdj(void)/释放动态数组内存if(m_pUnknownPoint!=NULL)delete m_pUnknownPoint;m_p

24、UnknownPoint=NULL;if(m_pKnownPoint!=NULL)delete m_pKnownPoint;m_pKnownPoint=NULL; if(m_pBearingObs!=NULL)delete m_pBearingObs;m_pBearingObs=NULL; if(m_pDistObs!=NULL)delete m_pDistObs;m_pDistObs=NULL;if (m_pZeroBearing!=NULL)delete m_pZeroBearing;m_pZeroBearing=NULL;/设置已知点个数和数组/注意：调用该函数后，原有数据被删除void

25、 CTraverseNetAdj:SetKnownPointSize(int size)if(m_pKnownPoint!=NULL)delete m_pKnownPoint;m_pKnownPoint=NULL;m_pKnownPoint=new CControlPointsize;/已知点数组m_iKnownPointCount=size;/已知点个数/设置未知点个数和数组/注意：调用该函数后，原有数据被删除void CTraverseNetAdj:SetUnknownPointSize(int size)if(m_pUnknownPoint!=NULL)delete m_pUnknown

26、Point;m_pUnknownPoint=NULL;m_iUnknownPointCount=size; m_pUnknownPoint=new CControlPointsize;/设置方向观测值个数和数组/注意：调用该函数后，原有数据被删除void CTraverseNetAdj:SetBearingObsSize(int size)if(m_pBearingObs!=NULL)delete m_pBearingObs;m_pBearingObs=NULL;m_pBearingObs=new CBearingObssize; /方向观测值数组m_iBearingObsCount=size

27、; /方向观测值个数/设置距离观测值个数和数组/注意：调用该函数后，原有数据被删除void CTraverseNetAdj:SetDistObsSize(int size)if(m_pDistObs!=NULL)delete m_pDistObs;m_pDistObs=NULL;m_pDistObs=new CDistObssize; /距离观测值数组m_iDistObsCount=size; /距离观测值个数/设置零方向观测值个数和数组/注意：调用该函数后，原有数据被删除void CTraverseNetAdj:SetZeroBearingSize(int size)if(m_pZeroBe

28、aring!=NULL)delete m_pZeroBearing;m_pZeroBearing=NULL;m_pZeroBearing=new CBearingObssize; /零方向观测值数组m_iZeroBearingCount=size; /零方向观测值个数/已知两个控制点，求P1-P2的方位角CAngle CTraverseNetAdj:Azi(const CControlPoint& P1, const CControlPoint& P2)CAngle angAzi; angAzi(RAD)=Azimuth(P1.X,P1.Y,P2.X,P2.Y);return angAzi;/

29、已知两个控制点，求P1-P2的距离double CTraverseNetAdj:HDist(const CControlPoint& P1, const CControlPoint& P2)return Dist(P1.X,P1.Y,P2.X,P2.Y);/字符串分割函数CString* CTraverseNetAdj:SplitString(CString str , char split, int &iSubStrs)int iPos = 0; /分割符位置 int iNums = 0; /分割符的总数 CString strTemp = str; CString strRight; /先

30、计算子字符串的数量 while (iPos != -1) iPos = strTemp.Find(split); if (iPos = -1) break; strRight = strTemp.Mid(iPos + 1, str.GetLength(); strTemp = strRight; iNums+; if (iNums = 0) /没有找到分割符 /子字符串数就是字符串本身 iSubStrs = 1; return NULL; /子字符串数组 iSubStrs = iNums + 1; /子串的数量= 分割符数量+ 1 CString* pStrSplit; pStrSplit =

31、 new CStringiSubStrs; strTemp = str; CString strLeft; for (int i = 0; i iNums; i+) iPos = strTemp.Find(split); /左子串 strLeft = strTemp.Left(iPos); /右子串 strRight = strTemp.Mid(iPos + 1, strTemp.GetLength(); strTemp = strRight; pStrSpliti = strLeft; pStrSplitiNums = strTemp; return pStrSplit;/导入观测数据函数b

32、ool CTraverseNetAdj:LoadObsData(const CString& strFileName)CStdioFile sf; /创建文件对象/以读的形式打开文件，如果打开失败则返回if(!sf.Open(strFileName, CFile:modeRead) return false;CString strLine;BOOL bEOF=sf.ReadString(strLine);/读取第一行，已知点个数 SetKnownPointSize(_ttoi(strLine);/根据已知点个数设置数组大小/开始读取已知点数据int n=0;/得到分割的子字符串数量/读取并保存

33、已知点数据 for(int i=0;im_iKnownPointCount;i+) sf.ReadString(strLine); CString *pstrData=SplitString(strLine,n); m_pKnownPointi.strName=pstrData0; m_pKnownPointi.strID=1; /将strID设为，即标记为坐标已知 m_pKnownPointi.X=_tstof(pstrData1); m_pKnownPointi.Y=_tstof(pstrData2); delete pstrData; pstrData=NULL;/开始读取未知点数据 s

34、f.ReadString(strLine);/未知点个数SetUnknownPointSize(_ttoi(strLine);/读取并保存未知点数据n=0;sf.ReadString(strLine);CString *pstrData=SplitString(strLine,n); for(int i=0;im_iUnknownPointCount;i+) m_pUnknownPointi.strName=pstrDatai; m_pUnknownPointi.strID=0; /将strID设为，即标记为坐标未知 delete pstrData; pstrData=NULL; /开始读取方

35、向观测值数据sf.ReadString(strLine);/方向观测值个数 SetBearingObsSize(_ttoi(strLine);/根据方向观测值个数设置数组大小n=0;for (int i=0;im_iBearingObsCount;i+)sf.ReadString(strLine);CString *pstrData=SplitString(strLine,n);m_pBearingObsi.cpStation=SearchPointUsingNAME(pstrData0);m_pBearingObsi.cpObject=SearchPointUsingNAME(pstrDat

36、a1);m_pBearingObsi.angleObsValue=_tstof(pstrData2);delete pstrData;pstrData=NULL;/根据方向观测值数组，将零方向观测值数据读入相应数组中int cnt=0; for(int i=0;im_iBearingObsCount;i+)if (m_pBearingObsi.angleObsValue(DEG)=0)m_iZeroBearingCount+;SetZeroBearingSize(m_iZeroBearingCount);for (int i=0;im_iBearingObsCount;i+)if (m_pBe

37、aringObsi.angleObsValue(DEG)=0)m_pZeroBearingcnt=m_pBearingObsi; cnt+;/开始读取距离观测值数据sf.ReadString(strLine);/距离观测值个数n=0;SetDistObsSize(_ttoi(strLine);/根据距离观测值个数设置数组大小 /读取并保存距离观测值数据 for(int i=0;im_iDistObsCount;i+) sf.ReadString(strLine); CString *pstrData=SplitString(strLine,n); m_pDistObsi.cpStart =Se

38、archPointUsingNAME(pstrData0); m_pDistObsi.cpEnd =SearchPointUsingNAME(pstrData1); m_pDistObsi.dDist =_tstof(pstrData2); delete pstrData; pstrData=NULL;sf.Close();return true;/根据点号从已知点数组中找到控制点，并返回该点的指针CControlPoint* CTraverseNetAdj:SearchKnownPointUsingNAME(CString strNAME)for(int i=0;im_iKnownPoint

39、Count;i+)if(strNAME=m_pKnownPointi.strName)return &m_pKnownPointi;return NULL;/根据点号从未知点数组中找到控制点，并返回该点的指针CControlPoint* CTraverseNetAdj:SearchUnknownPointUsingNAME(CString strNAME)for(int i=0;im_iUnknownPointCount;i+)if(strNAME=m_pUnknownPointi.strName)return &m_pUnknownPointi;return NULL;/根据点号从未知点和已

40、知点数组中找到控制点，并返回该点的指针CControlPoint* CTraverseNetAdj:SearchPointUsingNAME(CString strNAME)CControlPoint* pCP=NULL;pCP=SearchKnownPointUsingNAME(strNAME);if(pCP=NULL)pCP=SearchUnknownPointUsingNAME(strNAME);return pCP;/求近似坐标void CTraverseNetAdj:ApproximateXY()CAngle Ang;/两方向之间的夹角值for(int i=0;istrID=0&m_

41、pBearingObsi.cpStation-strID=1) /测站点已知，照准点未知for(int j=0;jstrName=(m_pBearingObsi.cpStation)-strName&m_pBearingObsj.cpObject-strID=1) /相同测站点，且照准点已知Ang(DEG)=m_pBearingObsi.angleObsValue(DEG)-m_pBearingObsj.angleObsValue(DEG); /夹角即为两方向值相减if(Ang(RAD)2*PI)EndBearing(RAD)=EndBearing(RAD)-2*PI;else if (End

42、Bearing(RAD)0)EndBearing(RAD)=EndBearing(RAD)+2*PI;for(int n=0;nstrName=(m_pBearingObsi.cpStation)-strName&(m_pDistObsn.cpEnd)-strName=(m_pBearingObsi.cpObject)-strName)(m_pBearingObsi.cpObject)-X=(m_pBearingObsi.cpStation)-X+m_pDistObsn.dDist*cos(EndBearing(RAD);(m_pBearingObsi.cpObject)-Y=(m_pBear

43、ingObsi.cpStation)-Y+m_pDistObsn.dDist*sin(EndBearing(RAD);(m_pBearingObsi.cpObject)-strID=1;break;/求未知改正数绝对值最大的double CTraverseNetAdj:MaxAbsDetX(const CMatrix& detX)double max=0;for(int i=0;imax)max=fabs(detX(i,0);return max;/把矩阵输出到文件中void CTraverseNetAdj:OutMatrixToFile(const CMatrix& mat,CStdioFi

44、le& SF) CString strLine,strTmp;for(int i=0;imat.Row();i+)strLine.Empty();for(int j=0;jmat.Col();j+)strTmp.Format(_T(%.4f ),mat(i,j);strLine=strLine+strTmp;SF.WriteString(strLine+_T(rn);/定权,SigmaAng为角度观测值先验精度（以秒为单位）,(方法仍需改进)/SigmaDist距离观测的精度（以厘米为单位）CMatrix CTraverseNetAdj:Weight(double SigmaAng, doub

45、le SigmaDist)int iTotalObsCount=m_iBearingObsCount+m_iDistObsCount; CMatrix P(iTotalObsCount, iTotalObsCount);P.Unit(); for(int i=0;im_iDistObsCount;i+)P(m_iBearingObsCount+i,m_iBearingObsCount+i)=(SigmaAng*SigmaAng)/(SigmaDist*SigmaDist*m_pDistObsi.dDist); return P;/计算定向零度方向角函数输入一个方向观察值返回测站点的零度定向角C

46、Angle CTraverseNetAdj:ZeroBearing(const CBearingObs m_BearingObs)CAngle alfa(0,RAD);/每个方向求出的坐标方位角CAngle Z(0,RAD);/零度定向角CAngle A(0,RAD);/每个方向观测值的照准角int j=0;for(int i=0;istrName=m_BearingObs.cpStation-strName)alfa=Azi(*m_pBearingObsi.cpStation,*m_pBearingObsi.cpObject);/求出测站点与照准点之间的近似坐标方位角A=alfa-m_pBe

47、aringObsi.angleObsValue;/每个方向求出来的坐标方位角减去方向观测值，得到每个方向观测值的照准角if(A(RAD)2*PI)A(RAD)=A(RAD)-2*PI;if(A(RAD)Y - BearingObs.cpStation-Y) / (Sjk*Sjk) / 1000;Bjk = -Ro * (BearingObs.cpObject-X - BearingObs.cpStation-X) / (Sjk*Sjk) / 1000; Aj = Ajk; Bj = Bjk; Ak = -Ajk; Bk = -Bjk;CAngle alfa=Azi(*BearingObs.cp

48、Station,*BearingObs.cpObject);/方向观测的方位角CAngle Z;/零方向的方位角CAngle LZ;/观测方向的方位角Z=ZeroBearing(BearingObs);LZ=BearingObs.angleObsValue+Z;if(LZ(RAD)2*PI)LZ(RAD)=LZ(RAD)-2*PI;if(LZ(RAD)X - DistObs.cpStart-X) / Sik; d = (DistObs.cpEnd-Y - DistObs.cpStart-Y) / Sik; Ci = -c; Di = -d; Ck = c; Dk = d; Li = (Dist

49、Obs.dDist - Sik) * 1000;/常数项以毫米为单位/搜索对应的零方向在零方向数组中的位置int CTraverseNetAdj:SearchZeroBearing(const CBearingObs &BearingObs)int j=0; for (int n=0;nstrName=BearingObs.cpStation-strName&m_pBearingObsn.angleObsValue(DEG)=0) for (j=0;jm_iZeroBearingCount;j+) if (m_pZeroBearingj.cpStation=m_pBearingObsn.cpS

50、tation) return j; return -1;/组成误差方程,B 为系数矩阵，L为常数项向量void CTraverseNetAdj:FormErrorEquations(CMatrix& B, CMatrix& L) int iObsCount;/观测值总个数int iUnknownCount;/误差方程中未知数的个数 iObsCount = m_iBearingObsCount + m_iDistObsCount; iUnknownCount = m_iUnknownPointCount*2+m_iZeroBearingCount; B.SetSize(iObsCount, iU

51、nknownCount);L.SetSize(iObsCount, 1); /计算方向观测值系数与常数项 double Aj, Bj, Ak, Bk, Li;for(int i=0;istrName)=NULL& SearchKnownPointUsingNAME(m_pBearingObsi.cpObject-strName)=NULL)/起点终点均为未知 for (j=0;jstrName=m_pUnknownPointj.strName) break; for (k=0;kstrName=m_pUnknownPointk.strName) break; B(i,j*2)=Aj; B(i,

52、j*2+1)=Bj; B(i,k*2)=Ak; B(i,k*2+1)=Bk; else if (SearchKnownPointUsingNAME(m_pBearingObsi.cpStation-strName)=NULL& SearchKnownPointUsingNAME(m_pBearingObsi.cpObject-strName)!=NULL)/起点为未知，终点为已知 for (j=0;jstrName=m_pUnknownPointj.strName) break; /*B(i,j)=Bi;*/ B(i,j*2)=Aj; B(i,j*2+1)=Bj; else if (Searc

53、hKnownPointUsingNAME(m_pBearingObsi.cpStation-strName)!=NULL& SearchKnownPointUsingNAME(m_pBearingObsi.cpObject-strName)=NULL)/对于起点为已知，终点为未知 for (k=0;kstrName=m_pUnknownPointk.strName) break; B(i,k*2)=Ak; B(i,k*2+1)=Bk; /计算边长观测值系数与常数项 double Ci, Di, Ck, Dk;for(int i=0;istrName)!=NULL&SearchKnownPoin

54、tUsingNAME(m_pDistObsi.cpEnd-strName)=NULL)/对于起点为已知，终点为未知 for (k=0;kstrName=m_pUnknownPointk.strName) break; B(i+m_iBearingObsCount,k*2)=Ck; B(i+m_iBearingObsCount,k*2+1)=Dk; else if (SearchKnownPointUsingNAME(m_pDistObsi.cpStart-strName)=NULL& SearchKnownPointUsingNAME(m_pDistObsi.cpEnd-strName)!=N

55、ULL) for (j=0;jstrName=m_pUnknownPointj.strName) break; /*B(i,j)=Bi;*/ B(i+m_iBearingObsCount,j*2)=Ci; B(i+m_iBearingObsCount,j*2+1)=Di; else if (SearchKnownPointUsingNAME(m_pDistObsi.cpStart-strName)=NULL& SearchKnownPointUsingNAME(m_pDistObsi.cpEnd-strName)=NULL) for (j=0;jstrName=m_pUnknownPointj

56、.strName) break; for (k=0;kstrName=m_pUnknownPointk.strName) break; B(i+m_iBearingObsCount,j*2)=Ci; B(i+m_iBearingObsCount,j*2+1)=Di; B(i+m_iBearingObsCount,k*2)=Ck; B(i+m_iBearingObsCount,k*2+1)=Dk; /严密平差（间接平差）计算主函数void CTraverseNetAdj:RigorousAdjust(const CString& strFileName)/开始输出间接平差的结果文件CStdioF

57、ile SF;CString strLine;setlocale(LC_ALL,); if(!SF.Open(strFileName, CFile:modeCreate|CFile:modeWrite) return; /开始写数据 SF.WriteString(_T(控制网间接平差结果n); /写已知点数据strLine.Format(_T(已知点个数：%dn),m_iKnownPointCount);SF.WriteString(strLine); for(int i=0;im_iKnownPointCount;i+) strLine.Format(_T(%s,%.4f,%.4fn), m

58、_pKnownPointi.strName, m_pKnownPointi.X,m_pKnownPointi.Y);SF.WriteString(strLine); /第一步：计算未知点近似坐标ApproximateXY(); /把未知点近似坐标写入文件strLine.Format(_T(未知点个数：%dn),m_iUnknownPointCount);SF.WriteString(strLine); for(int i=0;im_iUnknownPointCount;i+) strLine.Format(_T(%s,%.4f,%.4fn), m_pUnknownPointi.strName,

59、 m_pUnknownPointi.X,m_pUnknownPointi.Y);SF.WriteString(strLine); /第二步：计算方向观测值误差方程的系数和常数项 /第三步：计算边长观测值误差方程的系数和常数项 /第四步：组成误差方程/以上三步包含在FormErrorEquations函数中 /第五步：定权，得到权矩阵int iTotalObsCount=m_iBearingObsCount+m_iDistObsCount; CMatrix P(iTotalObsCount,iTotalObsCount); P=Weight(5, 0.5); /第六步：组成并解算法方程 CMat

60、rix B, L; CMatrix BT, N, InvN, detX, V;CMatrix W; /注意：detX，即坐标改正数的单位为mm(循环迭代发生改变的变量是什么,前面的系数计算有问题) do FormErrorEquations(B, L); BT = B; N = BT*P*B; InvN = N.Inv();W=BT*P*L; detX = InvN * (BT * P * L); for (int i=0;i 0.1); /第六步：精度评定 V = B * detX - L; /输出系数阵和常数项 SF.WriteString(_T(B矩阵：rn); OutMatrixToF