误差理论与测量平差上机指导书

1、误差理论与测量平差上机指导书辽宁工程技术大学测绘与地理科学学院测绘工程系误差理论与测量平差上机指导书目录实验1 矩阵加法与乘法运算1实验2 矩阵转置与求逆运算3实验3 误差椭圆元素计算11实验4 水准网间接平差程序设计12实验1 矩阵加法与乘法运算一、 实验名称：矩阵加法与乘法运算。二、 实验目的和任务：掌握矩阵加法与乘法通用程序的编写。三、 实验要求：1 每人独立编写出矩阵加法与乘法的程序，并上机调试通过；2 采用VC+6.0开发平台，C或者C+语言编写程序；3 写出矩阵运算的结果。四、 实验内容：1 矩阵加法的示例函数(C语言)void JZjiafa(double a1515,doubl

2、e b1515,double c1515,int m,int n)for (int i=0;i<=m-1;i+)for(int j=0;j<=n-1;j+)cij=aij+bij; return;2 矩阵乘法的示例程序(C语言)#include "stdafx.h"void matrixMultiply(double a1415,double b1513, double c1413,long m,long n,long k) for (long i = 0; i<= m-1; i+)for (long j=0; j<=k-1; j+)cij =0.0

3、;for (long q=0; q<=n-1;q+)cij = cij + aiq * bqj;return;int main(int argc, char* argv)long n,m,k,i,j;double a1415,c1413,b1513;FILE *stream;stream = fopen("矩阵输入.txt","r");fscanf(stream,"%ld %ld",&n,&m);for (i=0;i<n;i+)for(j=0;j<m;j+)fscanf(stream,"%l

4、f",&aij);fscanf(stream,"%ld %ld",&m,&k);for(i=0;i<m;i+)for(j=0;j<k;j+)fscanf(stream,"%lf",&bij);fclose(stream);matrixMultiply(a,b,c,4,5,3);stream = fopen("矩阵计算结果.txt","w");for (i=0;i<=3;i+)for(j=0;j<=2;j+)fprintf(stream,"%

5、16.7e ",cij);fprintf(stream,"n");fprintf(stream,"n");fclose(stream); return 0;实验2 矩阵转置与求逆运算一、 实验名称：矩阵转置与求逆运算。二、 实验目的和任务：掌握矩阵转置的编写，会调用矩阵求逆函数。三、 实验要求：1 每人独立编写出矩阵转置的程序，并上机调试通过；2 每人独立完成矩阵求逆函数的调用，并调试通过；3 采用VC+6.0开发平台，C或者C+语言编写程序；4 写出矩阵运算的结果。四、 实验内容：i. 矩阵的转置示例函数(C语言)double JZzhuan

6、zhi(double a1515, double b1515, int m,int n) for(int i=0;i<m;i+) for(int j=0;j<n;j+) bji=aij; return 0.0; ii. 矩阵求逆的示例函数(C语言)int invGJ(double *a,int n)int *is,*js,i,j,k,l,u,v;double d,p;is=(int *)malloc(n*sizeof(int);js=(int *)malloc(n*sizeof(int);for(k=0;k<=n-1;k+)d=0.0;for(i=k;i<=n-1;i+

7、)for(j=k;j<=n-1;j+)l=i*n+j;p=fabs(aij);if(p>d)d=p;isk=i;jsk=j; if(d+1.0=1.0)free(is);free(js);printf("error not invn");return (0);if(isk!=k)for(j=0;j<=n-1;j+)u=k*n+j;v=isk*n+j;p=akj;akj=aiskj;aiskj=p; if(jsk!=k) for(i=0;i<=n-1;i+) u=i*n+k;v=i*n+jsk; p=aik;aik=aijsk;aijsk=p; l=k

8、*n+k; akk=1.0/akk; for(j=0;j<=n-1;j+) if(j!=k) u=k*n+j;akj=akj*akk; for(i=0;i<=n-1;i+)if(i!=k) for(j=0;j<=n-1;j+) if(j!=k) u=i*n+j; aij=aij-aik*akj; for(i=0;i<=n-1;i+) if(i!=k) u=i*n+k;aik=-aik*akk; for(k=n-1;k>=0;k-) if(jsk!=k) for(j=0;j<=n-1;j+) u=k*n+j;v=jsk*n+j; p=akj;akj=ajskj

9、;ajskj=p; if(isk!=k)for(i=0;i<=n-1;i+)u=i*n+k;v=i*n+isk;p=aik;aik=aiisk;aiisk=p;free(is);free(js);return (1); int invGJ(double *a,int n)int *is,*js,i,j,k,l,u,v;double d,p;is=(int *)malloc(n*sizeof(int);js=(int *)malloc(n*sizeof(int);for(k=0;k<=n-1;k+)d=0.0;for(i=k;i<=n-1;i+)for(j=k;j<=n-

10、1;j+)l=i*n+j;p=fabs(aij);if(p>d)d=p;isk=i;jsk=j; if(d+1.0=1.0)free(is);free(js);printf("error not invn");return (0);if(isk!=k)for(j=0;j<=n-1;j+)u=k*n+j;v=isk*n+j;p=akj;akj=aiskj;aiskj=p; if(jsk!=k) for(i=0;i<=n-1;i+) u=i*n+k;v=i*n+jsk; p=aik;aik=aijsk;aijsk=p; l=k*n+k; akk=1.0/akk

11、; for(j=0;j<=n-1;j+) if(j!=k) u=k*n+j;akj=akj*akk; for(i=0;i<=n-1;i+)if(i!=k) for(j=0;j<=n-1;j+) if(j!=k) u=i*n+j; aij=aij-aik*akj; for(i=0;i<=n-1;i+) if(i!=k) u=i*n+k;aik=-aik*akk; for(k=n-1;k>=0;k-) if(jsk!=k) for(j=0;j<=n-1;j+) u=k*n+j;v=jsk*n+j; p=akj;akj=ajskj;ajskj=p; if(isk!

12、=k)for(i=0;i<=n-1;i+)u=i*n+k;v=i*n+isk;p=aik;aik=aiisk;aiisk=p;free(is);free(js);return (1);iii. 矩阵求逆函数的调用(C语言)#include <stdio.h>#include <stdlib.h>#include <math.h>int invGJ(double *a,int n);void main() int i,j;double *AA;/首先对二维指针Naa分配内存，采用C语言的方法/*AA=(double *)malloc(sizeof(dou

13、ble)*2);for(i=0;i<2;i+)AAi=(double *)mallo(sizeof(double)*2);*/首先对二维指针Naa分配内存，采用C+语言的方法AA=new double * 2;for(i=0;i<2;i+) AAi=new double2;double BB22=1,2,3,4;for(i=0;i<2;i+)for(j=0;j<2;j+)AAij=BBij;/调用矩阵求逆函数invGJ(AA,2);printf("矩阵AA的逆阵如下n");for(i=0;i<2;i+)for(j=0;j<2;j+)pri

14、ntf("%10.4lf",AAij);printf("n");double CC22;printf("AA与其逆阵的乘积如下(理论上是单位阵)n");for(i=0;i<2;i+)for(j=0;j<2;j+)CCij=0.0;for(int k=0;k<2;k+)CCij+=AAik*BBkj;printf("%10.4lf",CCij);printf("n");/C语言释放AA二维指针的方法/*for(i=0;i<2;i+)free(AAi);free(AA); *

15、/C+语言释放AA二维指针的方法for(i=0;i<2;i+)delete AAi;delete AA;实验3 误差椭圆元素计算一、 实验名称：误差椭圆元素计算。二、 实验目的和任务：掌握误差椭圆和相对误差椭圆元素的计算公式，并采用C或者C+语言变成实现。三、 实验要求：1 每人独立编写出误差椭圆和相对误差椭圆元素的计算程序，并调试通过；2 采用VC+6.0开发平台，C或者C+语言编写程序；3 写出计算的结果。4 本实验属于综合性，设计性实验，对学生的要求比较高，会综合使用矩阵加法，乘法以及转置和求逆的程序。5 注意事项：坐标方位角计算时必须考虑X与Y 所在的象限。四、 实验内容：参考教

16、材测量平差P103-P112。具体的例子可以参考P111页的例题，进行程序的调试。(一) 误差椭圆元素计算公式：(二) 相对误差椭圆元素计算公式：五、 例子（P111）在某三角网中插入P1及P2两个新点。设用间接平差法平差。平差之后这两点之间的协因数阵如下：根据以上的公式，分别计算出未知点P1和P2的误差椭圆元素以及这两点之间的相对误差椭圆元素。实验4 水准网间接平差程序设计一、 实验名称：水准网间接平差程序设计二、 实验目的：掌握用间接平差法对任意网形的水准网进行平差的算法设计以及程序编制；并学习采用读文件处理数据的方法。三、 实验任务：用C/C+编写水准网间接平差程序，并调试通过，用实测数

17、据实算分析；同时评定精度。四、 实验要求：1 个人独立编写程序，原始数据存放于文本文件或者数据库文件中；2 程序应该具有通用性，即任意网形都可以平差；3 用实测的数据进行计算分析。4 评定待定点的高程精度。五、 实验内容1 数据文件的编制格式总点数 未知点数 测段数（观测高差的个数）已知点1点名 高程已知点2点名 高程未知点3点名 0未知点4点名 0未知点5点名 0。起点点号 终点点号 观测高差 路线长度(km)起点点号 终点点号 观测高差 路线长度(km)起点点号 终点点号 观测高差 路线长度(km)。2 平差原理由观测值的起始和终点号或者水准网网形，形成误差方程的系数矩阵B(也叫设计矩阵)

18、，由观测的路线长度形成观测高差的权阵P(观测值独立，P为对角阵)，原理如下： (C为任意常数)3 法方程系数矩阵与闭合差的自动累加由于N和W具有可加性，即每读取一个高差观测值，即可得到一个误差方程的系数向量，然后累加到法方程系数矩阵与闭合差中；当高差观测值读取完毕的时候，法方程系数矩阵与闭合差也累加完毕；此时可用公式直接计算待定点高程的改正数，加上高程近似值，就得到了高程的平差值。4 精度评定 验后单位权中误差：待定点高程(未知参数)的协因数阵：待定点高程的方差：改正数V的协因数阵：高差平差值的协因数阵：5 例子1（P79）（1）in.txt文件4 3 51 A 237.4832 B 0.03

19、 C 0.04 D 0.01 2 5.835 3.52 3 3.782 2.71 3 9.640 4.04 3 7.384 3.01 4 2.270 2.5（2）源程序/Gckzwpc.h头文件代码class CGckzwpc public:CGckzwpc();virtual CGckzwpc();public:bool ReadData(CString filename); /读水准网平差数据文件void pc(); /平差函数void jdpd(); /精度评定函数int invGJ(double *a,int n); /求逆函数int nz,nw,ne,nn; /文件头信息：总点数，未

20、知点数，已知点数，测段数int n120,n220; /存放高差起点与终点的点号double H50,h50,W20,X20,B5020,V50,VPV;double *Nbb,S50,l50,P50,ph50;CString dm20; /控制点点名double SIG0,DX2020;/Gckzwpc.cpp代码bool CGckzwpc:ReadData(CString filename)int i;int MAXLINE =512;char buff513,ch115;CStdioFile fp;if( !fp.Open(filename,CFile:modeRead|CFile:ty

21、peText,NULL)AfxGetApp()->m_pMainWnd->MessageBox("数据文件不存在或数据文件错!","进程. . . . . .!",MB_OK|MB_ICONSTOP);_exit(1);return FALSE;fp.ReadString(buff,MAXLINE);sscanf(buff,"%d%d%d",&nz,&nw,&nn);ne=nz-nw;int dh;double gc;for(i=0;i<nz;i+)fp.ReadString(buff,MAX

22、LINE);sscanf(buff,"%d%s%lf",&dh,ch1,&gc);dmi=ch1;Hi=gc;for(i=0;i<nn;i+)fp.ReadString(buff,MAXLINE);sscanf(buff,"%d%d%lf%lf",&n1i,&n2i,&hi,&Si);fp.Close(); return TRUE;void CGckzwpc:pc()int i,j,k;/double Nbb12020;Nbb = new double * nw;for(i=0;i<nw;i+)

23、Nbbi = new double nw;for(i=0;i<nn;i+)Pi=0.0;li=0.0;for(j=0;j<nw;j+)Bij=0.0;for(i=0;i<nn;i+)if(n1i>ne) Bin1i-ne-1=-1;if(n2i>ne) Bin2i-ne-1=+1;Pi=1/Si;li=Hn1i-1+hi-Hn2i-1;for(i=0;i<nw;i+)for(j=0;j<nw;j+)Nbbij=0.0;for(k=0;k<nn;k+)Nbbij+=Bki*Bkj*Pk;for(i=0;i<nw;i+)Wi=0.0;for(k

24、=0;k<nn;k+)Wi+=Bki*Pk*lk;invGJ(Nbb,nw);for(i=0;i<nw;i+)Xi=0.0;for(k=0;k<nw;k+)Xi+=Nbbik*Wk;/精度评定函数void CGckzwpc:jdpd()int i,j,k;FILE *fp;fp=fopen("out.txt","w");VPV =0.0;for(i=0;i<nn;i+)Vi=0.0;for(k=0;k<nw;k+)Vi+=Bik*Xk;Vi+=-li;VPV+=Vi*Vi*Pi;phi=hi+Vi;/计算验后单位权中误差SI

25、G0=sqrt(VPV/(nn-nw);/计算未知参数的方差for(i=0;i<nw;i+)for(j=0;j<nw;j+)DXij=SIG0*SIG0*Nbbij;double QLL5050,SIGL5050;/计算观测值的平差值的中误差 QLL=B*inv(Nbb)*BT DLL=SIG0*SIG0*QLLfor(i=0;i<nn;i+)for(j=0;j<nn;j+)QLLij =0.0;SIGLij=0.0;for(k=0;k<nw;k+)for(int m=0;m<nw;m+)QLLij+=Nbbkm*Bik*Bjm;SIGLij=SIG0*sq

26、rt(QLLij);fprintf(fp,"*水准网间接平差结果*nn");fprintf(fp," 总点数%3d 未知点数%3d 测段数%3dn",nz,nw,nn);fprintf(fp,"n 验后单位权中误差%6.2lf(mm)nn",SIG0*1000);fprintf(fp,"n 起点号 终点号 观测高差(m) 路线长(km) 改正数(mm) 平差高差(m) 中误差(mm)n");for(i=0;i<nn;i+)fprintf(fp," %3d %3d %8.4lf %8.4lf %8.2

27、lf %8.4lf %8.2lfn",n1i,n2i,hi,Si,Vi*1000,phi,SIGLii*1000);fprintf(fp,"nn 已知点 已知高程(m)n");for(i=0;i<ne;i+)fprintf(fp,"%10s %10.4lf n",dmi,Hi);fprintf(fp,"n 待定点 平差高程(m) 高程中误差(mm)n");for(i=0;i<nw;i+)fprintf(fp," %10s %10.4lf %5.2lfn",dmi+ne,Xi,sqrt(DXii)*1000);for(i=0;i<nw;i+)delete Nbbi;delet