数据结构课程设计导线网平差程序的设计与实现

1、设计二：导线网平差程序的设计与实现一、设计目的立足于数据结构与测绘软件开发这一课程的课堂教学及其实验课程设计，为着实提高学生基于计算机辅助的方式切实解决工程实际问题的动手能力，通过本实习，一方面，使学生深入了解课堂所学知识，另一方面，通过实践掌握测绘行业软件设计与开发的基本方法，深刻掌握矩阵运算、曲线/曲面拟合的数值解法，掌握不同类型的典型测绘软件设计方法，使得学生初步具备编写测绘软件常用算法的能力以及开发中小规模测绘专业软件的能力。二、设计内容有导线网如图，观测了14条边长和16个转折角，已知测角精度，测边精度为。已知A、B、C、D、E、F点的坐标（无误差），如下表：表1 已知点数据点号X(

2、m)Y(m)A5256.9534520.068B5163.7524281.277C3659.3713621.210D4119.8793891.607E4581.1505345.292F4851.5545316.953表2 角度观测值编号角度观测值（ ）编号角度观测值（ ）1163 45 049169 10 30264 58 371098 22 043250 18 111194 53 504103 57 3412111 14 23583 08 051379 20 186258 54 1814268 06 047249 13 1715180 41 188207 32 3416103 23 08表3

3、 边长观测值编号边长观测值（m）编号边长观测值（m）1238.6198241.5602170.7599224.9963217.86910261.8264318.17311279.8405245.63512346.4436215.51413312.1097273.82914197.637试编程实现该导线网的平差，给出计算结果，并对其精度进行评价。三、关键问题描述3.1 未知点近似坐标计算平面控制网进行平差计算时需要计算未知点的近似坐标1.坐标计算公式1、2点的坐标已知，并观测了1-2、1-3的夹角，根据这些数据可以求出3号点坐标根据1、2两点的坐标，可以反算出1、2方向的方位角T12,3号点的坐

4、标为式子中S13为观测边长，为观测角度2.计算流程从读入的数据循环计算未知点的坐标，已计算出的坐标当做已知坐标的点处理参加下次计算，以此类推，逐步计算出未知点的坐标3.实现算法CMatrix CPlaneNetAdjust:XYJS()CMatrix _XYJS(Pnumber,2);double T12;for(int i=0;i0&xyk2.Y0) T12=GetT12(k1,k2);double s12=Gets12(k1,k2);double s13=Gets12(k1,k3);double T13=T12+guancejiaoi.Guancezhi;double dx=s13*cos

5、(T13);double dy=s13*sin(T13);xyk3.X=xyk1.X+dx; xyk3.Y=xyk1.Y+dy;for(int i=0;iPnumber;i+)double temp1=xyi.X;double temp2=xyi.Y;_XYJS.setValue(i,0,temp1);_XYJS.setValue(i,1,temp2);return _XYJS;3.2 误差方程列立1.理论分析平面控制网的误差方程都是非线性方程，必须引入参数近似值将误差方程线性化， 取的充分近似值 ，是微小量，在按台劳公式展开时可以略去二次和二次以上的项，而只取至一次项，于是可对非线性平差值观

6、测方程式线性化，于是有如下的式子对于观测角的改正数有对于边长观测值的改正数有2.实现算法如下：CMatrix CPlaneNetAdjust:B()CMatrix _B1(Lnumber,Pnumber*2);double a;double b;double c;double d;double m;double n;double m1;double n1;for(int i=0;iSnumber;i+)int k1=starti;int k2=endi;double dx=xyk2.X-xyk1.X;double dy=xyk2.Y-xyk1.Y;a=-dx/Gets12(k1,k2);b=-

7、dy/Gets12(k1,k2);if(k1knPnumber)_B1.setValue(i,2*k1,0);_B1.setValue(i,2*k1+1,0);else _B1.setValue(i,2*k1,a); _B1.setValue(i,2*k1+1,b);if(k2knPnumber)_B1.setValue(i,2*k2,0);_B1.setValue(i,2*k2+1,0);else_B1.setValue(i,2*k2,-a);_B1.setValue(i,2*k2+1,-b);for(int i=0;iTnumber;i+)const double p=206.265;in

8、t k1=cezhani;int k3=huoshii;int k2=qianshii;double dx12=xyk2.X-xyk1.X;double dy12=xyk2.Y-xyk1.Y;double dx13=xyk3.X-xyk1.X;double dy13=xyk3.Y-xyk1.Y;c=(p*dx13/Gets12(k1,k3)/Gets12(k1,k3)-p*dx12/Gets12(k1,k2)/Gets12(k1,k2);c=-c;d=-p*dy13/Gets12(k1,k3)/Gets12(k1,k3)+p*dy12/Gets12(k1,k2)/Gets12(k1,k2);d

9、=-d;m=-p*dy13/Gets12(k1,k3)/Gets12(k1,k3);m=-m;n=p*dx13/Gets12(k1,k3)/Gets12(k1,k3);n=-n;m1=p*dy12/Gets12(k1,k2)/Gets12(k1,k2);m1=-m1;n1=-p*dx12/Gets12(k1,k2)/Gets12(k1,k2);n1=-n1;if(k1=knPnumber)_B1.setValue(i+Snumber,2*k1,c);_B1.setValue(i+Snumber,2*k1+1,d);if(k2=knPnumber)_B1.setValue(i+Snumber,2

10、*k2,m1);_B1.setValue(i+Snumber,2*k2+1,n1);if(k3=knPnumber)_B1.setValue(i+Snumber,2*k3,m);_B1.setValue(i+Snumber,2*k3+1,n);CMatrix _B(Lnumber,2*(Pnumber-knPnumber);for(int i=0;i_B1.getRow();i+)for(int j=2*knPnumber;j2*Pnumber;j+)double temp=_B1.getValue(i,j);_B.setValue(i,(j-2*knPnumber),temp);return

11、 _B;3.3 法方程构建与解算1.理论分析误差方程系数构成法方程2.实现代码 计算LCMatrix CPlaneNetAdjust:L()CMatrix _L(Lnumber,1);double l;double s;double s0;for(int i=0;iSnumber;i+)int k1=starti;int k2=endi;s=Gets12(k1,k2);double dx=xyk2.X-xyk1.X;double dy=xyk2.Y-xyk1.Y;s0=sqrt(dx*dx+dy*dy);l=s-s0;_L.setValue(i,0,l*1000);/coutlendl;for

12、(int i=0;iTnumber;i+)const double pi=2*3.;double A12;double A13;double A;int k1=cezhani;int k2=huoshii;int k3=qianshii;A12=GetT12(k1,k2);A13=GetT12(k1,k3);A=GetA(k1,k2,k3);l=A13-A12;if(l=pi)l=l-pi;l=l-A;/coutrad_dms(A12) rad_dms(A13) rad_dms(A) rad_dms(l)endl;_L.setValue(i+Snumber,0,rad_dms(l);retur

13、n _L;计算权阵P CMatrix CPlaneNetAdjust:P()CMatrix _P(Lnumber,Lnumber);for(int i=0;iSnumber;i+)double temp=temp=Cjwucha*Cjwucha/(1*sqrt(bianchangi.Len)/(1*sqrt(bianchangi.Len); _P.setValue(i,i,temp);for(int i=Snumber;iLnumber;i+) _P.setValue(i,i,Cjwucha*Cjwucha/(Cjwucha*Cjwucha);return _P; 3.4 精度估计1.单位权中

14、误差间接平差与条件平差虽采用了不同的函数模型，但它们是在相同的最小乘原理下进行的，所以两法的平差结果总是相等的，这是因为在满足 条件下的V是唯一确定的，故平差值 不因方法不同而异。 单位权方差 的估值 计算式仍然是 除以其自由度，即2.协因数阵在间接平差中， 进一步有 3.误差椭圆参数点位误差椭圆参数E椭圆参数F4.实现代码void CPlaneNetAdjust:ErrorEllipse()CMatrix _Q=Q();const double PAI=3.;double m2=Cjwucha*Cjwucha;coutendlendl = 点位误差椭圆 =endl;cout 点名 点位误差

15、椭圆长半轴 椭圆短半轴 长轴方位角 endl;for(int i=0;i(Pnumber-knPnumber); i+)double mx2=_Q.getValue(2*i,2*i)*m2; /x坐标中误差的平方double my2=_Q.getValue(2*i+1,2*i+1)*m2; /y坐标中误差的平方double mxy=_Q.getValue(2*i,2*i+1)*m2; /xy坐标的协方差double K=sqrt(mx2-my2)*(mx2-my2)+4.0*mxy*mxy);double E=sqrt(0.5*(mx2+my2+K); /长轴double F=sqrt(0.5

16、*(mx2+my2-K); /短轴double A; /误差椭圆长轴的方位角if(fabs(mxy)my2) A=0.0;else A=0.5*PAI;else / mxy 0A=atan( (E*E-mx2)/mxy);if(A0.0)A+=PAI;coutendl;cout Pnamei+knPnumber Cjwucha*sqrt(_Q.getValue(2*i,2*i)+_Q.getValue(2*i+1,2*i+1) E F rad_dms(A);3.5 程序设计流程图近似坐标计算误差方程列立平差值计算精度评定四、算法运行结果与分析五、小结1.实习感受通过这次课程设计，我又对整本书有

17、了一个更深的理解。其实课程设计就是将我们所学的理论知识应用于实践的过程，在这一过程中，进一步掌握测量平差的基本原理和基本公式，并熟悉测量数据处理的基本技能和计算方法。这门课程与计算机程序完美地结合。这便要求我们在原有的解题思路中加入C+语言程序，并让它来帮助我们解决矩阵的复杂运算。既然用到了程序，我们就必须保证其运算的简洁性、正确性，尤其是在编写过程中要认真检查，为程序顺利运行打下基础。另外在各个子程序调用过程中，我们要充分考虑其顺序性并反复调试，以便得到理想结果。尽管在这次课程设计中遇到了很多困难，但我却得到了不少收获，并培养了自己正确应用公式、综合分析和解决问题的能力，同时也为今后步入社会

18、打下了一定的基础。另外，我们还要学会综合利用自身所学的知识，并将它们联系起来帮助自己有效地解决实际中的问题。总之，在这次课程设计中我不但过了比较充实的一周，还收获了不少知识。我或多或少明白，编程是一个循序渐进的过程，我想真的能够在这方面有所突破，就更要花时间多花工夫2.遇到的问题六、附录（完整代码）1.平面控制网平差类的设计在头文件PlaneNetAdjust.h下#include#includeMatrix.husing namespace std;struct XY int Index;double X;double Y;struct Bianchangstring startPoint;

19、 string endPoint; double Len;struct Guancejiaoint Index;string Cezhan;string Huozhan;string Qianzhan;double Guancezhi;class CPlaneNetAdjustpublic:CPlaneNetAdjust(void);CPlaneNetAdjust(void);void Getdata();/获取数据void Printdata();/数据显示void outdata();/数据输出CMatrix P();/计算权阵CMatrix L();/l的计算 CMatrix B();/

20、计算误差系数矩阵 CMatrix Q();/计算协因数阵CMatrix XYJS();/计算近似坐标CMatrix dx();/参数平差值CMatrix V();/残差计算void ErrorEllipse();/误差椭圆double Compute_T12(int k1,int k2);protected :double dms_rad(double a);double rad_dms(double a);int GetDianhao(string dianming);string GetDiaming(int dianhao);double Gets12(int k1,int k2);do

21、uble GetA(int k1,int k2,int k3);double GetT12(int k1,int k2);int Lnumber;int Snumber;int Tnumber;int Pnumber;int knPnumber;double Cjwucha;double Bcbiliwucha,Bcgudingwucha;XY *xy;Guancejiao *guancejiao;Bianchang *bianchang;CMatrix _XYJS;CMatrix _Q;CMatrix _B;CMatrix _P;CMatrix _L;CMatrix _dx;string *

22、Pname;int *start;int *end;int *cezhan;int *huoshi;int *qianshi;2.平面控制网类的实现在文件PlaneNetAdjust.cpp下#include stdafx.h#include PlaneNetAdjust.h#include#include#include#include using namespace std;CPlaneNetAdjust:CPlaneNetAdjust(void) Lnumber=0; Snumber=0; Tnumber=0; Pnumber=0; knPnumber=0;CPlaneNetAdjust

23、:CPlaneNetAdjust(void)delete Pname;delete xy;delete guancejiao;delete bianchang;delete start;delete end;double CPlaneNetAdjust:dms_rad(double a)/提取角度值的符号double sign=(a0.0) ? -1.0 : 1.0;a=fabs(a);/提取角度值的整度int d=(int)(a+0.00001)/10000.0);a=a-d*10000.0;if(a0.0) d=d-1; a=a+10000;/提取角度值的整分及秒值int m=(int)(

24、a+0.00001)/100.0);a=a-m*100;if(a0.0) m=m-1; a=a+100.0;a=sign*(d*3600.0+m*60.0+a)/.3;return a;/ 将角度的弧度值化为度分秒连写的角度（double 型） double CPlaneNetAdjust:rad_dms(double a)a=a*.3;double sign=(a0.0) ? -1.0 : 1.0;a=fabs(a);int d=(int)(a/3600.0+0.);a=a-d*3600.0;if(a0.0) d=d-1; a=a+3600.0; int m=(int)(a/60.0+0.0

25、001);a=a-m*60.0;if(a0.0) m=m-1; a=a+60.0; a=d*10000.0+m*100.0+a;return a*sign;void CPlaneNetAdjust:Getdata()string file=;ifstream inFile;coutendl;cout请输入文件的位置和文件名：file;coutendl;inFile.open(file,ios:in);if(!inFile)cout不能打开Lnumber;inFilePnumber;inFileknPnumber;inFileSnumber;inFileTnumber;inFile Cjwuch

26、aBcbiliwuchaBcgudingwucha;xy=new XYPnumber;Pname=new stringPnumber;for(int i=0;iPnamei;inFilexyi.X;inFilexyi.Y ;xyi.Index=i;bianchang=new BianchangSnumber;for(int i=0;ibianchangi.startPoint;inFilebianchangi.endPoint;inFilebianchangi.Len;start=new intSnumber;end=new intSnumber;for(int i=0;iSnumber;i+

27、)starti=GetDianhao(bianchangi.startPoint);endi=GetDianhao(bianchangi.endPoint);guancejiao=new GuancejiaoTnumber;double *dms=new doubleTnumber;for(int i=0;iguancejiaoi.Cezhan;inFileguancejiaoi.Huozhan;inFileguancejiaoi.Qianzhan;inFiledmsi;guancejiaoi.Guancezhi=dms_rad(dmsi);guancejiaoi.Index=i;cezhan

28、=new intTnumber;huoshi=new intTnumber;qianshi=new intTnumber;for(int i=0;iTnumber;i+)cezhani=GetDianhao(guancejiaoi.Cezhan);huoshii=GetDianhao(guancejiaoi.Huozhan);qianshii=GetDianhao(guancejiaoi.Qianzhan);cout 总点数：Pnumber 边长观测数：Snumber 角度观测数：Tnumberendl;/cout 测角中误差:Cjwucha 测边固定误差：Bcgudingwucha 测边比例

29、误差：Bcbiliwucha endlendl;coutendl;cout点的坐标（未知点坐标用0表示）：endl;cout点号 点名 X Yendl;coutsetprecision(3)setiosflags (ios:fixed);for(int i=0;iPnumber;i+)couti Pnamei xyi.X xyi.Y endl;cout边长观测：endl;cout起点 终点 观测值endl;for(int i=0;iSnumber;i+)coutbianchangi.startPoint bianchangi.endPoint bianchangi.Lenendl;cout角度

30、观测值：endl;cout测站 后视 前视 观测值endl;for(int i=0;iTnumber;i+)coutsetprecision(4)setiosflags (ios:fixed);coutguancejiaoi.Cezhan guancejiaoi.Huozhan guancejiaoi.Qianzhan dmsiendl;int CPlaneNetAdjust:GetDianhao(string dianming)for(int i=0;iPnumber;i+)if(Pnamei=dianming)return i; return -1;string CPlaneNetAdju

31、st:GetDiaming(int dianhao)return Pnamedianhao;double CPlaneNetAdjust:Gets12(int k1,int k2)for(int i=0;iSnumber;i+)if(k1=starti&k2=endi)return bianchangi.Len;if(k2=starti&k1=endi)return bianchangi.Len;return -1; /找不到符合条件的角，返回负值double CPlaneNetAdjust:GetA(int k1,int k2,int k3) for(int i=0;iTnumber;i+)

32、 if(k1=cezhani&k2=huoshii&k3=qianshii) return guancejiaoi.Guancezhi; if(k1=cezhani&k3=huoshii&k2=qianshii) return guancejiaoi.Guancezhi; return -1.0; /找不到符合条件的角，返回负值CMatrix CPlaneNetAdjust:P()CMatrix _P(Lnumber,Lnumber);for(int i=0;iSnumber;i+)double temp=temp=Cjwucha*Cjwucha/(1*sqrt(bianchangi.Len)

33、/(1*sqrt(bianchangi.Len); _P.setValue(i,i,temp);for(int i=Snumber;iLnumber;i+) _P.setValue(i,i,Cjwucha*Cjwucha/(Cjwucha*Cjwucha);return _P;double CPlaneNetAdjust:GetT12(int k1,int k2)const double Pi=3.;double dx=xyk2.X -xyk1.X ;double dy=xyk2.Y -xyk1.Y ;double T=atan2(dy,dx);if(T0)T=T+2.0*Pi;return

34、T;CMatrix CPlaneNetAdjust:XYJS()CMatrix _XYJS(Pnumber,2);double T12;for(int i=0;i0&xyk2.Y0) T12=GetT12(k1,k2);double s12=Gets12(k1,k2);double s13=Gets12(k1,k3);double T13=T12+guancejiaoi.Guancezhi;double dx=s13*cos(T13);double dy=s13*sin(T13);xyk3.X=xyk1.X+dx; xyk3.Y=xyk1.Y+dy;for(int i=0;iPnumber;i

35、+)double temp1=xyi.X;double temp2=xyi.Y;_XYJS.setValue(i,0,temp1);_XYJS.setValue(i,1,temp2);return _XYJS;CMatrix CPlaneNetAdjust:B()CMatrix _B1(Lnumber,Pnumber*2);double a;double b;double c;double d;double m;double n;double m1;double n1;for(int i=0;iSnumber;i+)int k1=starti;int k2=endi;double dx=xyk

36、2.X-xyk1.X;double dy=xyk2.Y-xyk1.Y;a=-dx/Gets12(k1,k2);b=-dy/Gets12(k1,k2);if(k1knPnumber)_B1.setValue(i,2*k1,0);_B1.setValue(i,2*k1+1,0);else _B1.setValue(i,2*k1,a); _B1.setValue(i,2*k1+1,b);if(k2knPnumber)_B1.setValue(i,2*k2,0);_B1.setValue(i,2*k2+1,0);else_B1.setValue(i,2*k2,-a);_B1.setValue(i,2*

37、k2+1,-b);for(int i=0;iTnumber;i+)const double p=206.265;int k1=cezhani;int k3=huoshii;int k2=qianshii;double dx12=xyk2.X-xyk1.X;double dy12=xyk2.Y-xyk1.Y;double dx13=xyk3.X-xyk1.X;double dy13=xyk3.Y-xyk1.Y;c=(p*dx13/Gets12(k1,k3)/Gets12(k1,k3)-p*dx12/Gets12(k1,k2)/Gets12(k1,k2);c=-c;d=-p*dy13/Gets12

38、(k1,k3)/Gets12(k1,k3)+p*dy12/Gets12(k1,k2)/Gets12(k1,k2);d=-d;m=-p*dy13/Gets12(k1,k3)/Gets12(k1,k3);m=-m;n=p*dx13/Gets12(k1,k3)/Gets12(k1,k3);n=-n;m1=p*dy12/Gets12(k1,k2)/Gets12(k1,k2);m1=-m1;n1=-p*dx12/Gets12(k1,k2)/Gets12(k1,k2);n1=-n1;if(k1=knPnumber)_B1.setValue(i+Snumber,2*k1,c);_B1.setValue(i+

39、Snumber,2*k1+1,d);if(k2=knPnumber)_B1.setValue(i+Snumber,2*k2,m1);_B1.setValue(i+Snumber,2*k2+1,n1);if(k3=knPnumber)_B1.setValue(i+Snumber,2*k3,m);_B1.setValue(i+Snumber,2*k3+1,n);CMatrix _B(Lnumber,2*(Pnumber-knPnumber);for(int i=0;i_B1.getRow();i+)for(int j=2*knPnumber;j2*Pnumber;j+)double temp=_B

40、1.getValue(i,j);_B.setValue(i,(j-2*knPnumber),temp);return _B;CMatrix CPlaneNetAdjust:L()CMatrix _L(Lnumber,1);double l;double s;double s0;for(int i=0;iSnumber;i+)int k1=starti;int k2=endi;s=Gets12(k1,k2);double dx=xyk2.X-xyk1.X;double dy=xyk2.Y-xyk1.Y;s0=sqrt(dx*dx+dy*dy);l=s-s0;_L.setValue(i,0,l*1

41、000);/coutlendl;for(int i=0;iTnumber;i+)const double pi=2*3.;double A12;double A13;double A;int k1=cezhani;int k2=huoshii;int k3=qianshii;A12=GetT12(k1,k2);A13=GetT12(k1,k3);A=GetA(k1,k2,k3);l=A13-A12;if(l=pi)l=l-pi;l=l-A;/coutrad_dms(A12) rad_dms(A13) rad_dms(A) rad_dms(l)endl;_L.setValue(i+Snumber,0,rad_dms(l);return _L;CMatrix CPlaneNetAdjust:Q()CMatrix N;N=B().transpose().operator*(P().operator*(B();return N.inverse();CMatrix CPlaneNetAdjust:dx()CMatrix W;W=B().transpose().operator*(P().operator*(L();_dx=Q().opera