误差理论与测量平差课程设计报告.docx
《误差理论与测量平差课程设计报告.docx》由会员分享，可在线阅读，更多相关《误差理论与测量平差课程设计报告.docx（30页珍藏版）》请在冰点文库上搜索。
误差理论与测量平差课程设计报告
n
一、目录----------------------------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-3
1.359
1.1
1
1-4
2.009
1.7
2
2-3
0.363
2.3
3
2-4
1.012
2.7
4
3-4
0.657
2.4
5
3-5
0.238
1.4
6
5-2
-0.595
2.6
7
（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.计算BTP
main（）
{
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，即NBB
main（）
{
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，即W
main（）
{
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）
#include"stdio.h"
#defineM3
voidmain（）
{
floatMAT[M][2*M];
floatMAT1[M][M];
floatt;
inti,j,k,l;
/***********************************************/
/*对矩阵进行初始化*/
for（i=0;ifor（j=0;j<2*M;j++）
MAT1[j]='\0';
/*对MAT1矩阵赋初值*/
for（i=0;ifor（j=0;jscanf（"%f",&MAT1[j]）;
/*打印目标矩阵?
*/
printf（"原矩阵为：
\n"）;
for（i=0;i{
for（j=0;jprintf（"%13.7f",MAT1[j]）;
printf（"\n"）;
}/********************************************/
/*对MAT1矩阵进行扩展,MAT1矩阵添加单位阵，由M*M变成2M*2M矩阵*/
for（i=0;ifor（j=0;j<2*M;j++）
if（jelseif（j==M+i）MAT[j]=1;
elseMAT[j]=0;
/*对M矩阵进行变换，使得前半部分矩阵成为单位阵，则*/
/*后半部分矩阵即为所求矩阵逆阵*/
for（i=0;i{
/*对第i行进行归一化*/
for（j=0;j<2*M;j++）
for（k=i+1;kMAT[j]=MAT[j]+MAT[k][j];
t=MAT;
for（j=i;j<2*M;j++）
MAT[j]=MAT[j]/t;
/*对矩阵进行行变换，使得第i列只有一个元素不为零，且为1*/
for（k=0;kif（k!
=i）
{
t=MAT[k];
for（l=i;l<2*M;l++）
MAT[k][l]=MAT[k][l]-MAT[l]*t;
}
}
/*将后半部分矩阵即所求矩阵逆阵存入MAT2矩阵。
*/
for（i=0;i{
for（j=0;jMAT1[j]=MAT[j+M];
printf（"\n"）;
}
/*********************************************/
/*输出所求的逆阵*/
printf（"逆阵为：
\n"）;
for（i=0;i{
for（j=0;jprintf（"%8.4f",MAT1[j]）;
printf（"\n"）;
}
}
4.求NBB-1W，即改正数x
main（）
{
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.计算Bx
main（）
{
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.计算VTP
main（）
{
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.计算VTPV
main（）
{
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及BT
B={{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.0}}
l={{0.0},
{0.0}，
{4.0}，
{3.0}，
{7.0}，
{2.0}，
{0.0}}
P={{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.3704,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.3846}}
BT={{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并运行程序得到BTP
BTP={{0.9091,0.0,0.4348,0.0,-0.4167,-0.7143,0.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.3846}}
2.替换第2个程序中的BTP、B并运行程序得到BTPB，即NBB
NBB={{2.4748，-0.4167，-0.7143},
{-0.4167，1.3753，0.0}，
{-0.7143，0.0，1.0989}}
3.替换第3个程序中的BTP、l并运行程序得到BTPl，即W
W={{-2.6063},
{4.0281}，
{1.4286}}
二、在矩阵的逆计算函数程序中进行以下操作
运行程序，按照提示及以上运算得到的矩阵NBB输入其元素，运行的结果即为NBB-1NBB-1={{0.5307，0.1608，0.3450},
{0.1608，0.7758，0.1045}，
{0.3450，0.1045，1.1342}}
三、再次在矩阵相乘计算函数的程序前提下，进行以下子程序的调用
1.替换第4个程序中的NBB-1、W并运行程序得到NBB-1W，即所选参数的改正数x
x={{-0.2426},
{2.8552}，
{1.1421}}
2.替换第5个程序中的B、x并运行程序得到Bx
Bx={{-0.2426},
{2.8552}，
{-0.2464}，
{2.8552}，
{3.0978}，
{1.3847}，
{-1.1421}}
3.根据V=Bx-l求出各观测值的改正数V，并写出VT，然后替换第6个程序中的VT、P并运行程序得到VTP
V={{-0.2426},
{2.8552}，
{-4.2426}，
{-0.1448}，
{-3.9022}，
{-0.6153}，
{-1.1421}}
VT=｛{-0.2426,2.8552,-4.2426,-0.1448,-3.9022,-0.6153,-1.1421}｝
VTP={{-0.2205,1.6794,-1.8447,-0.0536,1.6260,-0.4395,-0.4393}}
4.替换第7个程序中的VTP、V并运行程序得到VTPV
VTPV=19.7997
四、求出各个观测值平差值并按要求平定精度
X3=6.3748mX4=7.0279mX5=6.6122m
h1=1.3588mh2=2.0119mh3=0.3588mh4=1.0119mh5=0.6531mh6=0.2374mh7=-0.5961m
根据公式可求得单位权中误差为2.225mm
h34=X3-X4Q34=[1-10]NBB-1[1-10]T=0.9849
H3=X3Q34=[100]NBB-1[100]T=0.5307
H4=X4Q34=[010]NBB-1[010]T=0.7758
3、4点高差中误差为2.208mm
3号点高程中误差为1.621mm
4号点高程中误差为1.96mm
第七部分总结
通过这次误差理论与测量平差的课程设计，我又对整本书有了一个更深的理解。
其实课程设计就是将我们所学的理论知识应用于实践的过程，在这一过程中，进一步掌握测量平差的基本原理和基本公式，并熟悉测量数据处理的基本技能和计算方法。
或许我们已对《误差理论与测量平差》这本书的理论知识有了一定了解，但将它应用于实践依然是我们的一个难点，尤其是将这门课程与计算机程序完美地结合。
这便要求我们在原有的解题思路中加入C语言程序，并让它来帮助我们解决矩阵的复杂运算。
既然用到了程序，我们就必须保证其运算的简洁性、正确性，尤其是在编写过程中要认真检查，为程序顺利运行打下基础。
另外在各个子程序调用过程中，我们要充分考虑其顺序性并反复调试，以便得到理想结果。
尽管在这次课程设计中遇到了很多困难，但我却得到了不少收获，并培养了自己正确应用公式、综合分析和解决问题的能力，同时也为今后步入社会打下了一定的基础。
另外，我们还要学会综合利用自身所学的知识，并将它们联系起来帮助自己有效地解决实际中的问题。
总之，在这次课程设计中我不但过了比较充实的一周，还收获了不少知识。
#include
#include
#include
#include
#include
#include
usingnamespacestd;
classSZWPC
{
private:
intgcz_zs;//高差总数
intszd_zs;//总点数
intyz_szd_zs;//已知点数
doubles0;//单位权水准路线长度
doublem_pvv;//[pvv]
int*qsd_dh;//高差起点号
int*zd_dh;//高差终点号
char**dm;//点名地址数组
double*gcz;//观测值数组
double*szd_gc;//高程值数组
double*P;//观测值的权
double*BTPB,*BTPL;//法方程系数矩阵与自由项
double*dX;//高程改正数、平差值
double*V;//残差
doublezwc;//单位权中误差
public:
SZWPC（）;
~SZWPC（）;
intij（inti,intj）;//对称矩阵下标计算函数
boolinverse（doublea[],intn）;//对称正定矩阵求逆（仅存下三角元素）（参考他人）
voidinputdata（char*datafile）;//输入原始数据函数
intdm_dh（char*name）;//点名转点号
voidca_H0（）;//近似高程计算函数
voidca_BTPB（）;//法方程组成函数
voidca_dX（）;//高程平差值计算函数
voidprintresult（char*resultfile）;//精度估计与平差值输出函数
doubleca_V（）;//残差计算函数
voidzxecpc（char*resultfile）;//最小二乘平差函数
};
//////////构造函数
SZWPC:
:
SZWPC（）
{
gcz_zs=0;
szd_zs=0;
yz_szd_zs=0;
}
//////析构函数
SZWPC:
:
~SZWPC（）
{
if（gcz_zs>0）
{
delete[]qsd_dh;
delete[]zd_dh;
delete[]gcz;
delete[]P;
delete[]V;
}
if（szd_zs>0）
{
delete[]szd_gc;
delete[]BTPB;
delete[]BTPL;
delete[]dX;
for（inti=0;iif（dm[i]!
=NULL）delete[]（dm[i]）;
delete[]dm;
}
}
////////对称矩阵下标计算函数
intSZWPC:
:
ij（inti,intj）
{
return（i>=j）?
i*（i+1）/2+j:
j*（j+1）/2+i;
}
//对称正定矩阵求逆（仅存下三角元素）
boolSZWPC:
:
inverse（doublea[],intn）
{
double*a0=newdouble[n];
for（intk=0;k{
doublea00=a[0];
if（a00+1.0==1.0）
{
delete[]a0;
returnfalse;
}
for（inti=1;i{
doubleai0=a[i*（i+1）/2];
if（i<=n-k-1）a0[i]=-ai0/a00;
elsea0[i]=ai0/a00;
for（intj=1;j<=i;j++）
{
a[（i-1）*i/2+j-1]=a[i*（i+1）/2+j]+ai0*a0[j];
}
}
for（int
升级会员 