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

误差理论与测量平差课程设计报告课题：水准网严密平差及精度评定院校：系别：指导教师：班级：姓名：目录一、目录----------------------------1二、序言----------------------------2三、设计思路------------------------3四、程序流程图----------------------4五、程序及说明----------------------5六、计算结果-----------------------12七、总结---------------------------15第二部分序言1、课程设计的性质、目的和任务误差理论与测量平差是一门理论与实践并重的课程，其课程设计是测量数据处理理论学习的一个重要的实践环节，它是在我们学习了专业基础课“误差理论与测量平差基础”课程后进行的一门实践课程。其目的是增强我们对误差理论与测量平差基础理论的理解，牢固掌握测量平差的基本原理和基本公式，熟悉测量数据处理的基本技能和计算方法，灵活准确地应用于解决各类数据处理的实际问题，并能用所学的计算机理论知识，编制简单的计算程序。2、误差理论与测量平差课程和其它课程的联系和分工这次课程设计中所用的数学模型和计算方法是我们在误差理论与测量平差课程中所学的内容，所使用的C程序语言使我们在计算机基础课程中所学知识。误差理论与测量平差课程设计是测量平差和计算机程序设计等课程的综合实践与应用，同时也为我们今后步入工作岗位打下了一定基础。3、课程设计重点及内容本次课程设计重点是培养我们正确应用公式、综合分析和解决问题的能力，以及计算机编程能力。另外它要求我们完成1－2个综合性的结合生产实践的题目。如目前生产实践中经常用到的水准网严密平差及精度评定，边角网（导线）严密平差及精度评定等。此次我所选的课程设计课题是水准网严密平差及精度评定，其具体内容如下：根据题目要求，正确应用平差模型列出观测值条件方程、误差方程和法方程；解算法方程，得出平差后的平差值及各待定点的高程平差值；评定各平差值的精度和各高程平差值的精度。具体算例为：如图所示水准网，有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号点的高程中误差。第三部分设计思路一、解题步骤（1）此次设计我所采用的模型为间接平差模型，根据已知条件我们可知观测总数n=7，必要观测数t=3（则多余观测数r=n-t=4），因此我需先选定三个参数，即3、4、5点的最或然高程X3、X4、X5（X=X0+x，X30=6.375、X40=7.025、X50=6.611；其中X0为参数的近似值，x为其改正值）为参数。（2）列出条件方程，即将每一个观测量的平差值分别表达成所选参数的函数，H1+h1=X3、H1+h2=X4、H2+h3=X3、H2+h4=X4、X3+h5=X4、X3+h6=X5、X5+h7=H2；整理后得出误差方程，v1=x3、v2=x4、v3=x3-4、v4=x4-3、v5=-x3+x4-7、v6=-x3+x5-2、v7=-x5，即v=Bx-l的形式。（3）定权，令每千米的观测高差为单位权观测，即Pi=1/Si，从而可写出权阵P；根据误差方程式又可得其系数矩阵B和自由项l，并由它们组成法方程NBBx-W=0（其中NBB=BTPB，W=BTPl），法方程的个数等于所选参数的个数。（4）解算法方程，求出参数改正值x并计算参数的平差值X=X0+x。（5）由误差方程计算V，并求出观测量的平差值。为了检查平差计算的正确性，将所求的值代入条件方程，看其是否满足方程。（6）精度评定，计算单位权中误差，按照题设要求列出权函数式，再根据平差参数的协方差阵求出协因数，最后求出某段高差中误差，某些点的高程中误差。二、程序设计思想考虑到在解题过程中一些计算的复杂性，我们需借助一些技术将计算简单化，快捷化，因此在课程设计过程中，我们把一些C语言程序设计引入其中；通过一些简单、明了的程序及子函数调用，我们就可以很方便快捷的求出用笔算比较繁琐、费时的矩阵乘积、矩阵的逆（如BTPB、BTPl）等运算。第四部分程序流程图根据题目列出条件方程并写成误差方程的形式V=Bx-l↓确定权阵，根据误差方程得到矩阵B、l进而写出BT↓运用C程序语言求出BTP，进一步得到NBB=BTPB、W=BTPl并求出NBB-1↓用C程序求出参数的改正数x=NBB-1W↓根据C程序语言求Bx，进而由V=Bx-l写出各观测值的改正数↓根据L=L+V求出各观测值的平差值↓检验所求各值是否正确，若无误则往下进行，反之检查各步骤查出错误并改正↓由程序计算VTP进而求出VTPV，求单位权中误差，再根据权函数式、协因数传播定律评定各观测值及所求高程的精度第五部分程序及说明一、矩阵相乘计算函数#include“stdio.h”voidMatrix(a,b,m,n,k,c)intm,n,k;doublea[],b[],c[];{inti,j,l,u;for(i=0;i=m-1;i++)for(j=0;j=k-1;j++){u=i*k+j;c[u]=0.0;for(l=0;l=n-1;l++)c[u]=c[u]+a[i*n+l]*b[l*k+j];}return;}1.计算BTPmain(){inti,j;staticdoublea[3][7]=BT；staticdoublec[3][7],b[7][7]=P；Matrixmul(a,b,3,7,7,c);printf(“\n”);for(i=0;i=2;i++){for(j=0;j=6;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;}2.计算BTPB，即NBBmain(){inti,j;staticdoublea[3][7]=BTP；staticdoublec[3][3],b[7][3]=B；Matrixmul(a,b,3,7,3,c);printf(“\n”);for(i=0;i=2;i++){for(j=0;j=2;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;}3.计算BTPl，即Wmain(){inti,j;staticdoublea[3][7]=BTP；staticdoublec[3][1],b[7][1]=l；Matrixmul(a,b,3,7,1,c);printf(“\n”);for(i=0;i=2;i++){for(j=0;j=0;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;}二、矩阵的逆计算函数（求NBB-1）#includestdio.h#defineM3voidmain(){floatMAT[M][2*M];floatMAT1[M][M];floatt;inti,j,k,l;/***********************************************//*对矩阵进行初始化*/for(i=0;iM;i++)for(j=0;j2*M;j++)MAT1[j]='\0';/*对MAT1矩阵赋初值*/for(i=0;iM;i++)for(j=0;jM;j++)scanf(%f,&MAT1[j]);/*打印目标矩阵?*/printf(原矩阵为：\n);for(i=0;iM;i++){for(j=0;jM;j++)printf(%13.7f,MAT1[j]);printf(\n);}/********************************************//*对MAT1矩阵进行扩展,MAT1矩阵添加单位阵，由M*M变成2M*2M矩阵*/for(i=0;iM;i++)for(j=0;j2*M;j++)if(jM)MAT[j]=MAT1[j];elseif(j==M+i)MAT[j]=1;elseMAT[j]=0;/*对M矩阵进行变换，使得前半部分矩阵成为单位阵，则*//*后半部分矩阵即为所求矩阵逆阵*/for(i=0;iM;i++){/*对第i行进行归一化*/for(j=0;j2*M;j++)for(k=i+1;kM;k++)MAT[j]=MAT[j]+MAT[k][j];t=MAT;for(j=i;j2*M;j++)MAT[j]=MAT[j]/t;/*对矩阵进行行变换，使得第i列只有一个元素不为零，且为1*/for(k=0;kM;k++)if(k!=i){t=MAT[k];for(l=i;l2*M;l++)MAT[k][l]=MAT[k][l]-MAT[l]*t;}}/*将后半部分矩阵即所求矩阵逆阵存入MAT2矩阵。*/for(i=0;iM;i++){for(j=0;jM;j++)MAT1[j]=MAT[j+M];printf(\n);}/*********************************************//*输出所求的逆阵*/printf(逆阵为：\n);for(i=0;iM;i++){for(j=0;jM;j++)printf(%8.4f,MAT1[j]);printf(\n);}}4.求NBB-1W，即改正数xmain(){inti,j;staticdoublea[3][3]=NBB-1；staticdoublec[3][1],b[3][1]=W；Matrixmul(a,b,3,3,1,c);printf(“\n”);for(i=0;i=2;i++){for(j=0;j=0;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;}5.计算Bxmain(){inti,j;staticdoublea[7][3]=B；staticdoublec[7][1],b[3][1]=x；Matrixmul(a,b,7,3,1,c);printf(“\n”);for(i=0;i=6;i++){for(j=0;j=0;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;}6.计算VTPmain(){inti,j;staticdoublea[1][7]=VT；staticdoublec[1][7],b[7][7]=P；Matrixmul(a,b,1,7,7,c);printf(“\n”);for(i=0;i=0;i++){for(j=0;j=6;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;}7.计算VTPVmain(){inti,j;staticdoublea[1][7]=VTP；staticdoublec[1][1],b[7][1]=V；Matrixmul(a,b,1,7,1,c);printf(“\n”);for(i=0;i=0;i++){for(j=0;j=0;j++)printf(“%8.4f\t”,c[i][j];printf(“\n”);}printf(“\n”);return0;注：程序中有下划线部分在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.