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

一、目录 ----------------------------1

二、序言 ---------------------------- 2

三、设计思路 ------------------------ 3

四、程序流程图 ---------------------- 4

五、程序及说明 ---------------------- 5

六、计算结果 -----------------------12

七、总结 --------------------------- 15

第二部分序言

1、课程设计的性质、目的和任务

误差理论与测量平差是一门理论与实践并重的课程，其课程设计是测量数据处理理论学习的一个重要的实践环节，它是在我们学习了专业基础课“误差理论与测量平差基础”课程后进行的一门实践课程。其目的是增强我们对误差理论与测量平差基础理论的理解，牢固掌握测量平差的基本原理和基本公式，熟悉测量数据处理的基本技能和计算方法，灵活准确地应用于解决各类数据处理的实际问题，并能用所学的计算机理论知识，编制简单的计算程序。

2、误差理论与测量平差课程和其它课程的联系和分工

这次课程设计中所用的数学模型和计算方法是我们在误差理论与测量平差课程中所学的内容，所使用的C程序语言使我们在计算机基础课程中所学知识。误差理论与测量平差课程设计是测量平差和计算机程序设计等课程的综合实践与应用，同时也为我们今后步入工作岗位打下了一定基础。

3、课程设计重点及内容

本次课程设计重点是培养我们正确应用公式、综合分析和解决问题的能力，以及计算机编程能力。另外它要求我们完成1－2个综合性的结合生产实践的题目。如目前生产实践中经常用到的水准网严密平差及精度评定，边角网（导线）严密平差及精度评定等。此次我所选的课程设计课题是水准网严密平差及精度评定，其具体内容如下：

根据题目要求，正确应用平差模型列出观测值条件方程、误差方程和法方程；解算法方程，得出平差后的平差值及各待定点的高程平差值；评定各平差值的精度和各高程平差值的精度。

具体算例为：

如图所示水准网，有2个已知点，3个未知点，7个测段。各已知数据及观测值见下表（1）已知点高程H1= ， H2= （2）高差观测值(m)

（3）求各待定点的高程；3-4点的高差中误差；3号点、4号点的高程中误差。

第三部分设计思路

一、解题步骤

（1）此次设计我所采用的模型为间接平差模型，根据已知条件我们可知观测总数n=7，必要观测数t=3（则多余观测数r=n-t=4），因此我需先选定三个参数，即3、4、5点的最或然高程X3、X4、X5（X=X0+x，X30=、X40=、X50=；其中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）定权，令每千米的观测高差为单位权观测，即P i=1/S i，从而可写出权阵P；根据误差方程式又可得其系数矩阵B和自由项l，并由它们组成法方程N BB x-W=0（其中N BB=B T PB，W=B T Pl），法方程的个数等于所选参数的个数。

（4）解算法方程，求出参数改正值x并计算参数的平差值X=X0+x。

（5）由误差方程计算V，并求出观测量的平差值。为了检查平差计算的正确性，将所求的值代入条件方程，看其是否满足方程。

（6）精度评定，计算单位权中误差，按照题设要求列出权函数式，再根据平差参数的协方差阵求出协因数，最后求出某段高差中误差，某些点的高程中误差。

二、程序设计思想

考虑到在解题过程中一些计算的复杂性，我们需借助一些技术将计算简单化，快捷化，因此在课程设计过程中，我们把一些C语言程序设计引入其中；通过一些简单、明了的程序及子函数调用，我们就可以很方便快捷的求出用笔算比较繁琐、费时的矩阵乘积、矩阵的逆（如B T PB、B T Pl）等运算。

第四部分程序流程图

↓

第五部分程序及说明一、矩阵相乘计算函数

#include “”

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;c[u]=;

for(l=0;l<=n-1;l++)

c[u]=c[u]+a[i*n+l]*b[l*k+j];

}

return;

}

1.计算B T P

main()

{

int i,j;

static double a[3][7]=B T；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

}

2.计算B T PB，即N BB

main()

{

int i,j;

static double a[3][7]=B T P；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

}

3.计算B T Pl，即W

main()

{

int i,j;

static double a[3][7]=B T P；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

}

二、矩阵的逆计算函数（求N BB-1）

#include ""

#define M 3

void main()

{

float MAT[M][2*M];

float MAT1[M][M];

float t;

int i,j,k,l;

/***********************************************/

/*对矩阵进行初始化*/

for(i=0;i

for(j=0;j<2*M;j++)

MAT1[j]='\0';

/*对MAT1矩阵赋初值 */

for(i=0;i

for (j=0;j

scanf("%f",&MAT1[j]);

/*打印目标矩阵*/

printf("原矩阵为：\n");

for (i=0;i

{

for (j=0;j

printf("%",MAT1[j]);

printf("\n");

} /********************************************/

/*对MAT1矩阵进行扩展,MAT1矩阵添加单位阵，由M*M变成2M*2M矩阵 */

for(i=0;i

for(j=0;j<2*M;j++)

if (j

else if (j==M+i) MAT[j]=1;

else MAT[j]=0;

/*对M矩阵进行变换，使得前半部分矩阵成为单位阵，则 */

/*后半部分矩阵即为所求矩阵逆阵 */

for(i=0;i

{

/*对第i行进行归一化 */

for (j=0;j<2*M;j++)

for(k=i+1;k

MAT[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;k

if(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;j

MAT1[j]=MAT[j+M];

printf("\n");

}

/*********************************************/ /*输出所求的逆阵*/

printf("逆阵为：\n");

for(i=0;i

{

for(j=0;j

printf("%",MAT1[j]);

printf("\n");

}

4.求N BB-1W，即改正数x

main()

{

int i,j;

static double a[3][3]=N BB-1；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

}

5.计算Bx

main()

{

int i,j;

static double a[7][3]=B；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

}

6.计算V T P

main()

{

int i,j;

static double a[1][7]=V T；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

}

7.计算V T PV

main()

{

int i,j;

static double a[1][7]=V T P；

static double c[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(“%\t”,c[i][j];

printf(“\n”);

}

printf(“\n”);

return0;

注：程序中有下划线部分在C语言环境中运行时，需根据已知条件及所求结果进行替换！

第六部分计算结果根据条件方程及定权原则写出B、l、P及B T

B={{，，},

{，，}，

{，，}}

l={{},

{}，

{}}

P={{,,,,,,},

{,,,,,,}，

{,,,,,,}}

B T={{,,,,,,},

{,,,,,,}，

{,,,,,,}}

一、在矩阵相乘计算函数的程序前提下，进行以下子程序的调用

1.替换第1个程序中的B T 、P并运行程序得到B T P

B T P={{,,,,,,},

{,,,,,,}，

{,,,,,,}}

2.替换第2个程序中的B T P、B并运行程序得到B T PB，即N BB

N BB={{，，},

{，，}，

{，，}}

3.替换第3个程序中的B T P、l并运行程序得到B T Pl，即W

W={{},

{}，

{}}

二、在矩阵的逆计算函数程序中进行以下操作

运行程序，按照提示及以上运算得到的矩阵N BB输入其元素，运行的结果即为N BB-1 N BB-1={{，，}, {，，}，

{，，}}

三、再次在矩阵相乘计算函数的程序前提下，进行以下子程序的调用

1.替换第4个程序中的N BB-1、W并运行程序得到N BB-1W，即所选参数的改正数x

x={{},

{}，

{}}

2.替换第5个程序中的B、x并运行程序得到Bx

Bx={{},

{}，

{}}

3.根据V=Bx-l求出各观测值的改正数V，并写出V T，然后替换第6个程序中的V T、P并运行程序得到V T P V={{},

{}，

{}}

V T=｛{,,,,,,}｝

V T P={{,,,,,,}}

4.替换第7个程序中的V T P、V并运行程序得到V T PV

V T PV=

四、求出各个观测值平差值并按要求平定精度

X3= m X4= m X5= m

h1=h2=h3= h4= h5=h6= h7=

根据公式可求得单位权中误差为

h34=X3- X4Q34=[1 -1 0] N BB-1[1 -1 0]T=

H3=X3Q34=[1 0 0] N BB-1[1 0 0]T=

H4= X4 Q34=[0 1 0] N BB-1[0 1 0]T=

3、4点高差中误差为

3号点高程中误差为

4号点高程中误差为

第七部分总结

通过这次误差理论与测量平差的课程设计，我又对整本书有了一个更深的理解。其实课程设计就是将我们所学的理论知识应用于实践的过程，在这一过程中，进一步掌握测量平差的基本原理和基本公式，并熟悉测量数据处理的基本技能和计算方法。

或许我们已对《误差理论与测量平差》这本书的理论知识有了一定了解，但将它应用于实践依然是我们的一个难点，尤其是将这门课程与计算机程序完美地结合。这便要求我们在原有的解题思路中加入C语言程序，并让它来帮助我们解决矩阵的复杂运算。既然用到了程序，我们就必须保证其运算的简洁性、正确性，尤其是在编写过程中要认真检查，为程序顺利运行打下基础。另外在各个子程序调用过程中，我们要充分考虑其顺序性并反复调试，以便得到理想结果。

尽管在这次课程设计中遇到了很多困难，但我却得到了不少收获，并培养了自己正确应用公式、综合分析和解决问题的能力，同时也为今后步入社会打下了一定的基础。另外，我们还要学会综合利用自身所学的知识，并将它们联系起来帮助自己有效地解决实际中的问题。

总之，在这次课程设计中我不但过了比较充实的一周，还收获了不少知识。

#include

#include<>

#include

#include<>

using namespace std;

class SZWPC

{

private:

int gcz_zs; //高差总数

int szd_zs; //总点数

int yz_szd_zs; //已知点数

double s0;//单位权水准路线长度

double m_pvv; //[pvv]

int *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);//对称正定矩阵求逆(仅存下三角元素)(参考他人） 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_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(int i=0; i

if(dm[i]!=NULL)delete[](dm[i]);

delete []dm;

}

////// // 对称矩阵下标计算函数

int SZWPC::ij(int i,int j)

{

return (i>=j) i*(i+1)/2+j :j*(j+1)/2+i; }

// 对称正定矩阵求逆(仅存下三角元素)

bool SZWPC::inverse(double a[],int n)

double *a0=new double[n];

for(int k=0;k

{

double a00=a[0];

if(a00+==

{

delete []a0;

return false;

}

for(int i=1;i

{

double ai0 = a[i*(i+1)/2];

if(i<=n-k-1)a0[i]=-ai0/a00;

else a0[i]=ai0/a00;

for(int j=1;j<=i;j++)

{

a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+ai0*a0[j]; }

}

for(int i=1;i

{

a[(n-1)*n/2+i-1]=a0[i];

}

a[n*(n+1)/2-1]=a00;

}

delete []a0;

return true;

}

////////// 原始数据输入函数

void SZWPC::inputdata(char *datafile)

{

ifstream infile(datafile,ios::in);//声明输入句柄infile打开地址为datafile的文件并 if(!infile)

{

cerr<<" Open error!"<

}

infile>>gcz_zs>>szd_zs>>yz_szd_zs;

infile>>s0;

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 [gcz_zs];

V=new double [gcz_zs];

P=new double [gcz_zs];

dm=new char* [szd_zs];

for(int i=0;i

{

dm[i] = NULL;// dm_dh函数根据dm[i]是否为NULL确定dm[i]是否为点名地址

}

char buffer[128]; //临时数组，保存从文件中读到的点名

for(int i=0;i<=yz_szd_zs-1;i++)// 读取已知高程数据

{

infile>>buffer;

int c=dm_dh(buffer);

infile>>szd_gc[i];

}

for(int i=0;i

{

infile>>buffer; //读取高程起点名

qsd_dh[i]=dm_dh(buffer);

infile>>buffer;//读取高程终点

zd_dh[i]=dm_dh(buffer);

infile>>gcz[i]>>P[i]; //读取高差值与路线长度

P[i]=s0/P[i];//线路长转化为观测值的权

}

();

}

////////// 点名转点号，返回点名对应的点号

int SZWPC::dm_dh(char *name)

{

for(int i=0;i

{

if(dm[i]!=NULL)

{

if(strcmp(name,dm[i])==0)return i;//将待查点名与已经存入点名数组的点名比较，若存在返回点号

}

else

{

int len=strlen(name);//判断点名长度

dm[i]=new char[len+1];//为点名申请存储空间

strcpy(dm[i], name);//待查点名是一个新的点名，将新点名的地址放到dm数组中 return i;//返回点号

}

return -1; //dm数组已经存满，且没有待查点名

}

//////////高程近似值计算

void SZWPC::ca_H0()

{

for(int i=yz_szd_zs;i

for(int j=1;;j++)

{

int k=0; //计算出近似高程的点数

for(int i=0;i

{

int k1=qsd_dh[i]; //高差起点号

int k2=zd_dh[i]; //高差终点号

if(szd_gc[k1]> && szd_gc[k2]

{

szd_gc[k2]=szd_gc[k1]+gcz[i];//计算近似高程

k++;

}

else

if(szd_gc[k1]< && szd_gc[k2]>//k2点高程或高程近似值已知，k1点高程或高程近似值未知

{

szd_gc[k1]=szd_gc[k2]-gcz[i];//计算近似高程

k++;

}

if(k==(szd_zs-yz_szd_zs))break;//所有的近似高程计算完成，退出

}

//////////////////////////////

// 组成法方程

void SZWPC::ca_BTPB()

{

int t=szd_zs;

for(int i=0; i

for(int k=0; k

{

int i=qsd_dh[k];//获取点号

int j=zd_dh[k];//获取点号

double Pk=P[k];//获取权值

double lk=gcz[k]-(szd_gc[j]-szd_gc[i]);//获得第k个自由项 BTPL[i]-=Pk*lk;//获得法方程自由项

BTPL[j]+=Pk*lk;

BTPB[ij(i,i)]+=Pk;//获得法方程系数矩阵

BTPB[ij(j,j)]+=Pk;

BTPB[ij(i,j)]-=Pk;

}

/////////////////////////

// 高程平差值计算

void SZWPC::ca_dX()

{

for(int i=0;i

if(!inverse(BTPB,szd_zs))//矩阵求逆

{

cerr<<"法方程系数矩阵降秩!"<

}

for(int i=0; i

{

double xi=;

for(int j=0; j

{

xi+=BTPB[ij(i,j)]*BTPL[j];

}

dX[i]=xi;

szd_gc[i]+=xi;//计算高程平差值

}

/////////////////////

// 残差计算

double SZWPC::ca_V()

{

double pvv=;

for(int i=0;i<=gcz_zs-1;i++)

{

int k1=qsd_dh[i];

int k2=zd_dh[i];

V[i]=(szd_gc[k2]-szd_gc[k1]-gcz[i])*1000;

pvv+=V[i]*V[i]*P[i];

}

return(pvv);

}

////////////////////

// 原始数据和平差值输出

void SZWPC::printresult(char *resultfile)

{

double pvv=ca_V(); // 残差计算

ofstream outfile(resultfile,ios::out);//以输出方式打开文件，若文件不存在，创建文件 //

输出原始观测数据

outfile<

outfile<<" "<<"已知点数:"<

outfile<<"单位权水准路线长:"<

outfile<

for(int i=0;i<=yz_szd_zs-1;i++)

{

outfile<<" "<

outfile<

}

outfile<

outfile<<"起始点名"<<" "<<"终点点名"<<" "<<" 高差观测值(m)"<<" "<<" 两点间距离(km)"<

for(int i=0;i<=gcz_zs-1;i++)

{

outfile<<" "<

outfile<

}

zwc=sqrt(pvv/(gcz_zs-(szd_zs-yz_szd_zs)));//计算单位权中误差

outfile<