误差理论与测量平差课程设计报告

1、目录一、目录1二、序言2三、设计思路3四、程序流程图 4五、程序及说明5六、计算结果12七、总结15第二部分序言1、课程设计的性质、目的和任务误差理论与测量平差是一门理论与实践并重的课程，其课程设计是测量数据处理理论学习的一个重要的实践环节，它是在我们学习了专业基础课“误差理论与测量平差基础”课程后进行的一门实践课程。其目的是增强我们对误差理论与测量平差基础理论的理解，牢固掌握测量平差的基本原理和基本公式，熟悉测量数据处理的基本技能和计算方法，灵活准确地应用于解决各类数据处理的实际问题，并能用所学的计算机理论知识，编制简单的计算程序。2、误差理论与测量平差课程和其它课程的联系和分工这次课程设计

2、中所用的数学模型和计算方法是我们在误差理论与测量平差课程中所学的内容，所使用的C程序语言使我们在计算机基础课程中所学知识。误差理论与测量平差课程设计是测量平差和计算机程序设计等课程的综合实践与应用，同时也为我们今后步入工作岗位打下了一定基础。3、课程设计重点及内容本次课程设计重点是培养我们正确应用公式、综合分析和解决问题的能力，以及计算机编程能力。另外它要求我们完成12个综合性的结合生产实践的题目。如目前生产实践中经常用到的水准网严密平差及精度评定，边角网(导线)严密平差及精度评定等。此次我所选的课程设计课题是水准网严密平差及精度评定，其具体内容如下：根据题目要求，正确应用平差模型列出观测值条

3、件方程、误差方程和法方程；解算法方程，得出平差后 的平差值及各待定点的高程平差值；评定各平差值的精度和各高程平差值的精度。具体算例为：如图所示水准网，有 2个已知点，3个未知点，7个测段。各已知数据及观测值见下表(1)已知点高程 H1=5.016m , H2=6.016m(2)高差观测值(m)端点号高差观测值m测段距离km序号1-31.3591.111-42.0091.722-30.3632.332-41.0122.743-40.6572.453-50.2381.465-2-0.5952.67(3)求各待定点的高程；3-4点的高差中误差；3号点、4号点的高程中误差。第三部分设计思路一、解题步骤

4、(1)此次设计我所采用的模型为间接平差模型，根据已知条件我们可知观测总数n=7,必要观测数t=3 (则多余观测数r=n-t=4 ),因此我需先选定三个参数， 即3、4、5点的最或然高程 头3、X、X5 (X=X0+x, X=6.375、 X0=7.025、X0=6.611 ;其中尤为参数的近似值，x为其改正值)为参数。(2)列出条件方程，即将每一个观测量的平差值分别表达成所选参数的函数，H+hi=X” H+h2=X、H+h3=X3H2+h4=Xl、X3+h5=X4X3+h6=Xs、+h7=H；整理后得出误差方程，Vl=X3、V2=X4、V3=X3-4、V4=X4-3、V5=X3+X4-7 V

5、6=-X 3+X5-2 V7=-X 5 ,即 v=Bx-1 的形5。(3)定权，令每千米的观测高差为单位权观测，即Pi=1/Si,从而可写出权阵 P;根据误差方程式又可得其系数矩阵B和自由项l ,并由它们组成法方程 NBbx-W=0 (其中NBb=BtPB, W=BPl),法方程的个数等于所选参 数的个数。(4)解算法方程，求出参数改正值x并计算参数的平差值 X=X，+Xo(5)由误差方程计算 V,并求出观测量的平差值。为了检查平差计算的正确性，将所求的值代入条件方程， 看其是否满足方程。(6)精度评定，计算单位权中误差，按照题设要求列出权函数式，再根据平差参数的协方差阵求出协因数， 最后求出

6、某段高差中误差，某些点的高程中误差。二、程序设计思想考虑到在解题过程中一些计算的复杂性，我们需借助一些技术将计算简单化，快捷化，因此在课程设 计过程中，我们把一些 C语言程序设计引入其中；通过一些简单、明了的程序及子函数调用，我们就可以很 方便快捷的求出用笔算比较繁琐、费时的矩阵乘积、矩阵的逆(如BTPR BTP1)等运算。第四部分程序流程图根据题目列出条件方程并写成误差方程的形式V=Bx-l确定权阵，根据误差方程得到矩阵B l进而写出BT运用C程序语言求出 BTP,进一步得到 Nbb=BfPB W=BPl并求出NBb-1用C程序求出参数的改正数 x=Nbb-W根据C程序语言求Bx,进而由V=

7、Bx-l写出各观测值的改正数根据L=L+V求出各观测值的平差值检验所求各值是否正确，若无误则往下进行，反之检查各步骤查出错误并改正由程序计算VTP进而求出VPV,求单位权中误差，再根据权函数式、协因数传播定律评定各观测值及所求高程的精度第五部分程序及说明一、矩阵相乘计算函数#includestdio.h void Matrix(a,b,m,n,k,c)int m,n,k;double a,b口,c；int i,j,l,u;for(i=0;i=m-1;i+)for(j=0;j=k-1;j+)u=i*k+j;cu=0.0;for(l=0;l;Matrixmul(a,b,3,7,7,c); prin

8、tf(n );for(i=0;i=2;i+) for(j=0;j=6;j+)printf(%8.4ft ，cij;printf(”)； printf(n );return。； 2 .计算BTPB,即Nh main() int i,j；static double a37=BTf；static double c33,b73=B;Matrixmul(a,b,3,7,3,c); printf(n );for(i=0;i=2;i+) for(j=0;j=2;j+)printf(%8.4ft ，cij;printf(”)； printf(n );return。； 3 .计算BTPl ,即W main()

9、int i,j；static double a37= 旺；static double c31,b71=l_；Matrixmul(a,b,3,7,1,c);printf(n );for(i=0;i=2;i+) for(j=0;j=0;j+)printf(%8.4ft ，cij;printf(”)；printf( “ n ” );return0;二、矩阵的逆计算函数(求NBB-1 )#include stdio.h #define M 3 void main()float MATM2*M;float MAT1MM;float t;int i,j,k,l;/*/* 对矩阵进行初始化*/ for(i=

10、0;iM;i+)for(j=0;j2*M;j+)MAT1j=0;/*对MAT1矩阵赋初值*/for(i=0;iM;i+)for (j=0;jM;j+)scanf(%f,&MAT1j);/* 打印目标矩阵?*/printf( 原矩阵为：n);for (i=0;iM;i+)for (j=0;jM;j+) printf(%13.7f,MAT1j); /*printf(n);*/ /*对MAT1矩阵进行扩展，MAT1矩阵添加单位阵，由 M*M变成2M*2M矩阵*/for(i=0;iM;i+)for(j=0;j2*M;j+)if (jM) MATj=MAT1j;else if (j=M+i) MATj=

11、1;else MATj=0;/*对M矩阵进行变换，使得前半部分矩阵成为单位阵，则 */*后半部分矩阵即为所求矩阵逆阵*/for(i=0;iM;i+)/*对第i行进行归一化*/for (j=0;j2*M;j+)for(k=i+1;kM;k+)MATj=MATj+MATkj;t=MAT;for(j=i;j2*M;j+)MATj=MATj/t;/*对矩阵进行行变换，使得第 i列只有一个元素不为零，且为 1*/for(k=0;kM;k+)if(k!=i)t=MATk;for (l=i;l2*M;l+)MATkl=MATkl-MATl*t;/*将后半部分矩阵即所求矩阵逆阵存入MAT2巨阵。*/for(i

12、=0;iM;i+)for(j=0;jM;j+)MAT1j=MATj+M;printf(n);/*/*输出所求的逆阵*/printf(逆阵为：n);for(i=0;iM;i+)for(j=0;jM;j+)printf(%8.4f,MAT1j);printf(n);4 .求NBb-1 W即改正数xmain() int i,j； static double a33=NBb-1 ;static double c31,b31=WyMatrixmul(a,b,3,3,1,c); printf(n );for(i=0;i=2;i+) for(j=0;j=0;j+)printf(%8.4ft ，cij;pri

13、ntf(n ); printf(n );return。； 5 .计算Bxmain() int i,j； static double a73=B；static double c71,b31=x；Matrixmul(a,b,7,3,1,c); printf(n );for(i=0;i=6;i+) for(j=0;j=0;j+) printf(%8.4ft ，cij;printf(”)； printf( 行”)； return。； 6.计算VTP main() int i,j； static double a17= V; static double c17,b77=P;Matrixmul(a,b,1

14、,7,7,c); printf(n );for(i=0;i=0;i+) for(j=0;j=6;j+) printf(%8.4ft ，cij;printf(”)； printf(n );return。； 7.计算VTPV main() int i,j； static double a17=VTP；static double c11,b71=V；Matrixmul(a,b,1,7,1,c); printf(n );for(i=0;i=0;i+) for(j=0;j=0;j+) printf(%8.4ft ，cij;printf(”)； printf(n );return。； 注：程序中有下划线部

15、分在C语言环境中运行时，需根据已知条件及所求结果进行替换！第六部分计算结果根据条件方程及定权原则写出B、 l 、 P 及BTB=1.0，0.0 ，0.0,0.0，1.0 ，0.0，1.0，0.0 ，0.0，0.0，1.0 ，0.0，-1.0，1.0， 0.0 ，-1.0，0.0， 1.0 ，0.0，0.0 ，-1.0l=0.0,0.0，4.0，3.0，7.0，2.0，0.0P=0.9091,0.0,0.0,0.0,0.0,0.0,0.0, 0.0,0.5882,0.0,0.0,0.0,0.0,0.0，0.0,0.0,0.4348,0.0,0.0,0.0,0.0，0.0,0.0,0.0,0.37

16、04,0.0,0.0,0.0，0.0,0.0,0.0,0.0,0.4167,0.0,0.0，0.0,0.0,0.0,0.0,0.0,0.7143,0.0，0.0,0.0,0.0,0.0,0.0,0.0,0.3846BT=1.0,0.0,1.0,0.0,-1.0,-1.0,0.0, 0.0,1.0,0.0,1.0,1.0,0.0,0.0，0.0,0.0,0.0,0.0,0.0,1.0,-1.0一、在矩阵相乘计算函数的程序前提下，进行以下子程序的调用1. 替换第 1 个程序中的BT 、 P 并运行程序得到BT PBTP=0.9091,0.0,0.4348,0.0,-0.4167,-0.7143,0

17、.0, 0.0,0.5882,0.0,0.374,0.4167,0.0,0.0，0.0,0.0,0.0,0.0,0.0,0.7143,-0.38462. 替换第2个程序中的BTP、B并运行程序得到BTPB,即NbbNBB=2.4748 ， -0.4167 ， -0.7143,-0.4167 ， 1.3753 ， 0.0 ，-0.7143 ， 0.0 ， 1.09893. 替换第 3 个程序中的BTP、 l 并运行程序得到BTPl ，即WW=-2.6063,4.0281 ，1.4286二、在矩阵的逆计算函数程序中进行以下操作运行程序，按照提示及以上运算得到的矩阵NBb输入其元素，运行的结果即为N

18、bb-1 N3b-1=0.5307 ,0.1608 ,0.3450,0.1608 ， 0.7758， 0.1045 ，0.3450 ， 0.1045， 1.1342三、再次在矩阵相乘计算函数的程序前提下，进行以下子程序的调用1 .替换第4个程序中的Nbb1、W并运行程序得到 NBB1W即所选参数的改正数Xx=-0.2426,2.8552 ，1.14212 . 替换第 5 个程序中的B、 X 并运行程序得到BXBX=-0.2426,2.8552，-0.2464，2.8552，3.0978，1.3847，-1.14213 .根据V=Bx-l求出各观测值的改正数V并写出V,然后替换第6个程序中的 V

19、 P并运行程序得到 VTPV=-0.2426,2.8552，-4.2426，-0.1448，-3.9022，-0.6153，-1.1421VT= -0.2426,2.8552,-4.2426,-0.1448,-3.9022,-0.6153,-1.1421VTP=-0.2205,1.6794,-1.8447,-0.0536,1.6260,-0.4395,-0.43934 .替换第7个程序中的VTP、V并运行程序得到 VTPVVTPV=19.7997四、求出各个观测值平差值并按要求平定精度X3=6.3748 mX4=7.0279 mX5=6.6122 mh1=1.3588m h2=2.0119m

20、h3=0.3588m h4=1.0119m h5=0.6531m h6=0.2374m h7=-0.5961m根据公式可求得单位权中误差为2.225mmh34= X3- X 4 Q34=1 -1 0 NBB-1 1 -1 0T=0.9849H3=X3 Q34=1 0 0 N BB-1 1 0 0 T =0.5307H4= X4 Q 34=0 1 0 N BB-10 1 0 T=0.77583、 4 点高差中误差为2.208mm3 号点高程中误差为1.621mm4 号点高程中误差为1.96mm第七部分总结通过这次误差理论与测量平差的课程设计，我又对整本书有了一个更深的理解。其实课程设计就是将我们

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

22、自己正确应用公式、综合分析和解决问题的能力，同时也为今后步入社会打下了一定的基础。另外，我们还要学会综合利用自身所学的知识，并将它们联系起来帮助自己有效地解决实际中的问题。总之，在这次课程设计中我不但过了比较充实的一周，还收获了不少知识。#include#include#include#include#include#includeusing namespace std;class SZWPCprivate:int gcz_zs; /高差总数int szd_zs; /总点数int yz_szd_zs; /已知点数double s0;/单位权水准路线长度高差起点号高差终点号点名地址数组观测值数组

23、高程值数组观测值的权double m_pvv; /pvvint *qsd_dh; /int *zd_dh; /char *dm; /double *gcz; /double *szd_gc; /double *P; /double *BTPB,*BTPL; / 法方程系数矩阵与自由项double *dX; /高程改正数、平差值double *V; / 残差double zwc; /单位权中误差public:SZWPC();SZWPC();int ij(int i,int j);/对称矩阵下标计算函数bool inverse(double a,int n);/对称正定矩阵求逆( 仅存下三角元素)

24、( 参考他人)void inputdata(char *datafile);/输入原始数据函数int dm_dh(char *name); /点名转点号void ca_H0(); /近似高程计算函数void ca_BTPB(); / 法方程组成函数void ca_dX(); /高程平差值计算函数void printresult(char *resultfile); /精度估计与平差值输出函数double ca_V(); /残差计算函数void zxecpc(char *resultfile);/最小二乘平差函数;/ / 构造函数SZWPC:SZWPC()gcz_zs=0;szd_zs=0;yz

25、_szd_zs=0;/ / 析构函数SZWPC:SZWPC()if(gcz_zs0)delete qsd_dh;delete zd_dh;delete gcz;delete P;delete V;if(szd_zs0)delete szd_gc;delete BTPB;delete BTPL;delete dX;for(int i=0; i=j)? i*(i+1)/2+j :j*(j+1)/2+i;/ 对称正定矩阵求逆（ 仅存下三角元素）bool SZWPC:inverse(double a,int n) double *a0=new doublen;for(int k=0;kn;k+)dou

26、ble a00=a0;if(a00+1.0=1.0)delete a0;return false;for(int i=1;in;i+)double ai0 = ai*(i+1)/2;if(i=n-k-1)a0i=-ai0/a00;else a0i=ai0/a00;for(int j=1;j=i;j+)a(i-1)*i/2+j-1=ai*(i+1)/2+j+ai0*a0j;for(int i=1;in;i+)a(n-1)*n/2+i-1=a0i;an*(n+1)/2-1=1.0/a00;delete a0;return true;原始数据输入函数/ void SZWPC:inputdata(ch

27、ar *datafile)ifstream infile(datafile,ios:in);/声明输入句柄infile 打开地址为datafile 的文件并if(!infile)cerr Open error!gcz_zsszd_zsyz_szd_zs;infiles0;szd_gc=new double szd_zs;dX=new double szd_zs;BTPB=new double szd_zs*(szd_zs+1)/2;BTPL=new double szd_zs;qsd_dh=new int gcz_zs;zd_dh=new int gcz_zs;gcz=new double g

28、cz_zs;V=new double gcz_zs;P=new double gcz_zs;dm=new char* szd_zs;for(int i=0;iszd_zs;i+)dmi = NULL;/ dm_dh函数根据dmi是否为NULL确定dmi是否为点名地址char buffer128; / 临时数组，保存从文件中读到的点名for(int i=0;ibuffer;int c=dm_dh(buffer);infileszd_gci;for(int i=0;ibuffer; /qsd_dhi=dm_dh(buffer);infilebuffer;/zd_dhi=dm_dh(buffer);

29、infilegcziPi; /Pi=s0/Pi;/读取高程终点读取观测数据读取高程起点名读取高差值与路线长度线路长转化为观测值的权infile.close();/ 点名转点号，返回点名对应的点号int SZWPC:dm_dh(char *name)for(int i=0;iszd_zs;i+)if(dmi!=NULL)if(strcmp(name,dmi)=0)return i;/将待查点名与已经存入点名数组的点名比较，若存在返回点号判断点名长度为点名申请存储空间elseint len=strlen(name);/dmi=new charlen+1;/strcpy(dmi, name);/待查

30、点名是一个新的点名，将新点名的地址放到dm数组中return i;/返回点号return -1; /dm数组已经存满，且没有待查点名/ 高程近似值计算void SZWPC:ca_H0()for(int i=yz_szd_zs;iszd_zs;i+)szd_gci=-10000.9;/为计算机设置辨别未知高程点的标志for(int j=1;j+)int k=0; /计算出近似高程的点数for(int i=0;i-10000.0 & szd_gck2-10000.0)/k1点高程或高程近似值已知，k2 点高程或高程近似值未知szd_gck2=szd_gck1+gczi;/计算近似高程k+;else

31、if(szd_gck1-10000.0)/k2点高程或高程近似值已知， k1 点高程或高程近似值未知szd_gck1=szd_gck2-gczi;/计算近似高程k+;所有的近似高程计算完成，退出 if(k=(szd_zs-yz_szd_zs)break;/ 组成法方程void SZWPC:ca_BTPB()int t=szd_zs;for(int i=0; iszd_zs*(szd_zs+1)/2; i+) BTPBi=0.0;/赋初值for(int i=0; iszd_zs; i+) BTPLi=0.0;/赋初值for(int k=0; kgcz_zs; k+)int i=qsd_dhk;/

32、获取点号int j=zd_dhk;/获取点号double Pk=Pk;/获取权值double lk=gczk-(szd_gcj-szd_gci);/获得第 k 个自由项BTPLi-=Pk*lk;/获得法方程自由项BTPLj+=Pk*lk;BTPBij(i,i)+=Pk;/获得法方程系数矩阵BTPBij(j,j)+=Pk;BTPBij(i,j)-=Pk;/ 高程平差值计算void SZWPC:ca_dX()for(int i=0;iyz_szd_zs;i+) BTPBij(i,i)=1.0e30;/处理已知点if(!inverse(BTPB,szd_zs)/矩阵求逆退出程序cerr 法方程系数矩

33、阵降秩!endl;/ 矩阵为奇异矩阵，无法求逆exit(0);/for(int i=0; iszd_zs; i+)/计算高程改正数double xi=0.0;for(int j=0; jszd_zs; j+)xi+=BTPBij(i,j)*BTPLj;dXi=xi;szd_gci+=xi;/计算高程平差值/ 残差计算double SZWPC:ca_V()double pvv=0.0;for(int i=0;i=gcz_zs-1;i+)int k1=qsd_dhi;int k2=zd_dhi;Vi=(szd_gck2-szd_gck1-gczi)*1000;pvv+=Vi*Vi*Pi;retur

34、n(pvv);/ 原始数据和平差值输出void SZWPC:printresult(char *resultfile)double pvv=ca_V(); / 残差计算ofstream outfile(resultfile,ios:out);/以输出方式打开文件，若文件不存在，创建文件/输出原始观测数据outfileendl 观测总数：gcz_zs 总点数：szd_zs;outfile 已知点数：yz_szd_zsendlendl;outfile单位权水准路线长：s0(km)endl;outfileendl=已知高程=endl; 输出原始观测数据已知点点号、高程for(int i=0;i=yz

35、_szd_zs-1;i+)outfile dmi;outfilesetiosflags(ios：fixed);outfilesetw(10)setprecision(4)szd_gciendl;outfileendlendl=高 差 观 测 值=endlendl;/输出原始观测数据高程观测值与路线长outfile起始点名 终点点名 高差观测值(m) 两点间距离(km)endl;for(int i=0;i=gcz_zs-1;i+)outfile dmqsd_dhisetw(9)dmzd_dhi;outfilesetiosflags(ios：fixed);outfilesetw(16)setpre

36、cision(4)gczi;outfilesetiosflags(ios：fixed);outfilesetw(16)setprecision(4)s0/Piendl;zwc=sqrt(pvv/(gcz_zs-(szd_zs-yz_szd_zs);/计算单位权中误差outfileendl= 单位权中误差=endl;/ 输出 单位权中误差outfileendl0=zwc(mm)endl;outfileendl=高 程 平 差 值 及 其 精 度=endlendl;/ 输出高程平差值及其精度outfile 点名近似高程(m) 改正数(mm)高程平差值(m) 中误差(mm)endl;for(int i=0; iszd_zs; i+)(outfilesetw(2)dmi;double dx=dXi;double qii=BTPBij(i,i);outfilesetiosflags(ios:fixed);/用定点格式显示浮点数outfilesetw(12)setprecision(4)szd_gci-dx;outfilesetiosflags(ios:fixed);outfilesetw(15)setprecision