水准网程序设计.docx
- 文档编号:25549507
- 上传时间:2023-06-09
- 格式:DOCX
- 页数:23
- 大小:158.31KB
水准网程序设计.docx
《水准网程序设计.docx》由会员分享,可在线阅读,更多相关《水准网程序设计.docx(23页珍藏版)》请在冰豆网上搜索。
水准网程序设计
误差理论与测量平差基础
(MATLAB)
实习报告
学号:
姓名:
班级:
专业:
测绘工程
课程名称:
误差理论与测量平差基础
任课老师:
20XX年5月
一、任务概述
利用MATLAB或者C++编程间接平差程序,通过该程序读取观测数据文件,并计算出平差结果。
C++
//水准网平差.cpp:
定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include
#include"math.h"
#include
#include"string.h"
usingnamespacestd;
boolinverse(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(inti=1;i { a[(n-1)*n/2+i-1]=a0[i]; } a[n*(n+1)/2-1]=1.0/a00; } delete[]a0; returntrue; } intmatrix_s(inti,intj) { return(i>=j)? i*(i+1)/2+j: j*(j+1)/2+i; } classCLevelingAdjust { public: CLevelingAdjust(); virtual~CLevelingAdjust(); intm_Lnumber;//高差总数 intm_Pnumber;//总点数 intm_knPnumber;//已知点数 intunPnumber;//未知点个数 doublem_pvv;//[pvv] FILE*resultfp;//文件指针,输出计算结果 doublem_Sigma;//验前单位权中误差(粗差探测、闭合差计算用) int*StartP;//高差起点号 int*EndP;//高差终点号 char**Pname;//点名地址数组 double*L;//观测值数组 double*Height;//高程值数组 double*P;//观测值的权 double*ATPA,*ATPL;//法方程系数矩阵与自由项 double*dX;//参数平差值(高程改正数) double*V;//残差数组 doublem_mu;//验后单位权中误差 voidInputdata(char*datafile);//输入原始数据函数 voidPrintdata();//输出原始数据函数 intGetStationNumber(char*name);//获取点号函数 voidca_H0();//近似高程计算函数 voidca_ATPA();//法方程组成函数 voidca_dX();//平差值计算函数 voidPrintResult();//精度估计与平差值输出函数 doubleca_V();//残差计算函数 }; CLevelingAdjust: : CLevelingAdjust() { m_Lnumber=0; m_Pnumber=0; } CLevelingAdjust: : ~CLevelingAdjust() { if(m_Lnumber>0) { delete[]StartP; delete[]EndP; delete[]L; delete[]P; delete[]V; } if(m_Pnumber>0) { delete[]Height; delete[]ATPA; delete[]ATPL; delete[]dX; for(inti=0;i if(Pname[i]! =NULL)delete[](Pname[i]); delete[]Pname; } } //原始数据输入函数 voidCLevelingAdjust: : Inputdata(char*datafile) { FILE*fr; if((fr=fopen(datafile,"r"))==NULL) { cout<<"\n数据文件打不开! "< exit(0); } fscanf(fr,"%d%d%d",&m_Lnumber,&m_Pnumber,&m_knPnumber); unPnumber=m_Pnumber-m_knPnumber; Height=newdouble[m_Pnumber]; dX=newdouble[m_Pnumber]; ATPA=newdouble[m_Pnumber*(m_Pnumber+1)/2]; ATPL=newdouble[m_Pnumber]; StartP=newint[m_Lnumber]; EndP=newint[m_Lnumber]; L=newdouble[m_Lnumber]; V=newdouble[m_Lnumber]; P=newdouble[m_Lnumber]; //fscanf(fp,"%lf",&m_Sigma); Pname=newchar*[m_Pnumber]; for(inti=0;i { //GetStationNumber函数根据Pname[i]是否为NULL //确定Pname[i]是否为点名地址 Pname[i]=NULL; } charbuffer[100]; //读取已知高程数据 for(inti=0;i<=m_knPnumber-1;i++) { fscanf(fr,"%s",buffer); intc=GetStationNumber(buffer); fscanf(fr,"%lf",&Height[c]); } //读取观测数据 for(inti=0;i { fscanf(fr,"%s",buffer);//读取高程起点名 StartP[i]=GetStationNumber(buffer); if(StartP[i]<0) { fprintf(resultfp,"\n数据文件出错: "); fprintf(resultfp,"\n第%d个高差的起始点名为\"%s\"",i+1,buffer); fclose(resultfp); exit(0); } fscanf(fr,"%s",buffer);//读取高程终点 EndP[i]=GetStationNumber(buffer); if(EndP[i]<0) { fprintf(resultfp,"\n数据文件出错: "); fprintf(resultfp,"\n第%d个高差终点的点名为\"%s\"",i+1,buffer); fclose(resultfp); exit(0); } fscanf(fr,"%lf%lf",&L[i],&P[i]);//读取高差值与路线长度 P[i]=1.0/P[i]; } fclose(fr); } //原始数据写到结果文件 voidCLevelingAdjust: : Printdata() { inti; fprintf(resultfp,"\n观测值总数: %d总点数: %d已知点数: %d\n", m_Lnumber,m_Pnumber,m_knPnumber); //fprintf(resultfp,"\n验前单位权中误差: %lf",m_Sigma); fprintf(resultfp,"\n\n已知高程: \n"); for(i=0;i<=m_knPnumber-1;i++) { fprintf(resultfp,"\n%5d%8s",i+1,Pname[i]); fprintf(resultfp,"%10.4lf",Height[i]); } fprintf(resultfp,"\n\n高差观测值: \n"); for(i=0;i<=m_Lnumber-1;i++) { fprintf(resultfp,"\n%5d%8s%8s",i+1,Pname[StartP[i]],Pname[EndP[i]]); fprintf(resultfp,"%12.4lf%10.3lf",L[i],1.0/P[i]); } } //点名存贮,返回点名对应的点号 intCLevelingAdjust: : GetStationNumber(char*name) { for(inti=0;i { if(Pname[i]! =NULL) { //将待查点名与已经存入点名数组的点名比较 if(strcmp(name,Pname[i])==0)returni; } else { //待查点名是一个新的点名,将新点名的地址放到Pname数组中 intlen=strlen(name); Pname[i]=newchar[len+1]; strcpy(Pname[i],name); returni; } } return-1;//Pname数组已经存满,且没有待查点名 } //高程近似值计算 voidCLevelingAdjust: : ca_H0() { for(inti=m_knPnumber;i intjj=0;//计算出近似高程的点数 for(intii=1;;ii++) { for(inti=0;i { intk1=StartP[i];//高差起点号 intk2=EndP[i];//高差终点号 if(Height[k1]>-9999.0&&Height[k2]<-9999.0) { Height[k2]=Height[k1]+L[i]; jj++; } elseif(Height[k1]<-9999.0&&Height[k2]>-9999.0) { Height[k1]=Height[k2]-L[i]; jj++; } } if(jj==(m_Pnumber-m_knPnumber))break; if(ii>(m_Pnumber-m_knPnumber)) { fprintf(resultfp,"\n\n下列点无法计算出概略高程: \n"); for(inti=0;i { if(Height[i]<-9999.0)printf("\n%s",Pname[i]); } cout<<"近似高程计算失败! "< fclose(resultfp); exit(0); } } } voidCLevelingAdjust: : ca_ATPA() { intt=m_Pnumber; for(inti=0;i for(inti=0;i for(inti=0;i { ATPA[matrix_s(i,i)]=1.0E08; } for(intk=0;k { inti=StartP[k]; intj=EndP[k]; doublePk=P[k]; doubleLk=L[k]-(Height[j]-Height[i]); ATPL[i]-=Pk*Lk; ATPL[j]+=Pk*Lk; ATPA[matrix_s(i,i)]+=Pk; ATPA[matrix_s(j,j)]+=Pk; ATPA[matrix_s(i,j)]-=Pk; } } voidCLevelingAdjust: : ca_dX() { if(! inverse(ATPA,m_Pnumber)) { cout<<"法方程系数矩阵降秩! "< fclose(resultfp); exit(0); } for(inti=0;i { doublexi=0.0; for(intj=0;j { xi+=ATPA[matrix_s(i,j)]*ATPL[j]; } dX[i]=xi; Height[i]+=xi; } } doubleCLevelingAdjust: : ca_V() { doublepvv=0.0; for(inti=0;i<=m_Lnumber-1;i++) { intk1=StartP[i]; intk2=EndP[i]; V[i]=Height[k2]-Height[k1]-L[i]; pvv+=V[i]*V[i]*P[i]; } m_mu=sqrt(pvv/(unPnumber)); return(m_mu); } voidCLevelingAdjust: : PrintResult() { fprintf(resultfp,"\n\n====高程平差值及其精度====\n"); fprintf(resultfp,"\n点名近似高程改正数高程平差值中误差\n"); for(inti=0;i { fprintf(resultfp,"\n%5s",Pname[i]); doubledx=dX[i]; doubleqii=ATPA[matrix_s(i,i)]; fprintf(resultfp,"%12.4lf%9.4lf%12.4lf%9.4lf", Height[i]-dx,dx,Height[i],sqrt(qii)*m_mu); } } int_tmain(intargc,_TCHAR*argv[]) { char*datafile="算例\\data.txt"; char*resultfile="算例\\result.txt"; printf("\n\n水准网平差\n"); printf("数据文件: %s\n",datafile); printf("结果文件: %s\n",resultfile); FILE*fp=fopen(resultfile,"w"); if(fp==NULL) { cout<<"创建结果文件失败! "< return-1; } CLevelingAdjustlevel; level.resultfp=fp; level.Inputdata(datafile); level.Printdata(); level.ca_H0(); level.ca_ATPA(); level.ca_dX(); level.ca_V(); level.PrintResult(); fclose(fp); cout<<"解算成功! "< charjj; cin>>jj; return0; MATLAB 图一 图二 图三 图四 三、水准网图 四、输入的数据格式 数据格式为TXT文件,如图所示: TXT文件格式说明: (1)第一行格式 第一行分别表示观测个数5个,水准点数4个,未知点3个,已知点1个,所有数据用英文逗号隔开 (2)已知点数据格式 第二行开始是已知点点号和高程,一行列一个已知点点号和高程,由于该水准网只有一个已知点,所有只能列出一行。 图中表示已知点点号为1,高程为237.483m (3)测站起始点号格式 该部分表示测站的起始点点号 (4)测站终点点号格式 该部分表示测站的终点点号 (5)高差格式 该部分表示各测站的高差 (6)距离格式 该部分表示各测站的距离S 五、流程图 六、附件代码 functionSDJianJiePingCha() [FileName,PathName]=uigetfile('*.txt','打开水准观测数据');%打开文件 f=csvread(strcat(PathName,FileName));%打开文件并存在矩阵f中 point=f(1,2);%获取所有水准点个数 n=f(1,1);%获得观测个数n t=f(1,3);%获得必要观测个数t y=f(1,4);%获得已知点个数y XX=zeros(point,1);%初始化XX阵等于0,方便下面把已知点高程和未知点参数估值放到XX阵 B=zeros(n,t);%初始化B阵,方便下面求V=Bx-l中的系数阵B; forj=1: y XX(j,1)=f(j+1,2);%把已知点高程放到XX阵中 end data=f((2+y): end,: );%从文件中获取观测数据,并放到data阵中 h=data(: 3);%从data中获取观测高差,并放到h阵中 P=zeros(n);%初始化权阵P forj=1: n P(j,j)=10/data(j,4);%以10km观测值为单位权误差计算权阵P end fori=1: n%通过循环求B阵 point1=data(i,1);%获取某个测站的起始点号 point2=data(i,2);%获取某个测站的终点点号 ifpoint1>y&&point2>y%当某测站起始点和终点高程都未知时,求B阵第i行 B(i,point1-y)=-1; B(i,point2-y)=1; elseifpoint1<=y&&point2>y%当起始点高程已知和终点高程未知时,求B阵第i行 B(i,point2-y)=1; XX(point2,1)=XX(point1,1)+h(i,1);%求第i个参数估值 elseifpoint1>y&&point2<=y%当起始点高程未知和终点高程已知时,求B阵第i行 B(i,point1-y)=-1; XX(point1,1)=XX(point2,1)-h(i,1);%求第i个参数估值 end end l=zeros(n,1);%初始化小l阵,方便下面求V=Bx-l中的系数阵l; fori=1: n%通过循环求小l point1=data(i,1); point2=data(i,2); l(i,1)=-(XX(point2,1)-XX(point1,1)-h(i,1)); end %带入间接平差数学模型公式进行计算: r=n-t;%求多余观测数 N=B'*P*B; W=B'*P*l; x=N\W; X=XX((y+1): end,1)+x; V=B*x-l; L=h+V; a0=sqrt(V'*P*V/r); Qxx=inv(N); Dxx=a0*a0*inv(N); %输出计算结果: disp('参数改正数: ') x=x' disp('参数平差值: ') X=X' disp('观测值改正数: ') V=V' disp('观测值平差值: ') L=L' disp('协方差阵: ') Dxx disp('单位权方差: ') a0 disp('协因数阵: ') Qxx B l P N W end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 水准 程序设计